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ABSTRACT 


A model of the behavior of a packed bed nuclear reactor fuel element is developed. It 
is capable of predicting the temperature, pressure and velocity fields as a function of position 
throughout the fuel element in both transient and steady state conditions. It is the starting 
point for the design of a real time analysis module for a reactor power controller. 


The fuel element consists of a packed bed of fuel particles between two concentrically 
mounted retention elements. It is cooled by hydrogen flowing radially inward through the 
bed. 


The model is based on the fundamental principles of mass, momentum and energy con- 
servation. The balances are applied to a two dimensional array of control volumes, using the 
lumped parameter approach, to generate sets of simultaneous linear semi-implicit finite 
difference equations. The Pressure Implicit with Splitting of Operators (PISO) algorithm is 
then used to advance the model variables in time. An energy balance derived from a single 
node model of a fuel particle is performed on the solid phase. 


The model code is applied to a series of steady state and transient problems, varying 
the peak power density from 0 to 2.1 GW/m'. The model predicts significant axial flow in 
the fuel particle bed. If the proper flow distribution is obtained along the inlet retention ele- 
ment, it is not maintained in the particle bed. This redistribution leads to a 300 K tempera- 
ture variation along the outlet plenum at high power. The magnitude of the redistribution is a 
function of power level. A design change in the shape of the inlet frit may give a more 
uniform temperature distribution at the full power condition. The model also predicts signifi- 
cant effects from conduction and radiation in the solid phase. 
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CHAPTER ] 
INTRODUCTION 


1.1 PROBLEM STATEMENT 


MIT, in conjunction with the Sandia National Laboratory, is developing an automated 
control system for a space based reactor. The proposed reactor consists of several annular 
packed bed fuel elements cooled by hydrogen flowing radially inward through the fuel par- 
ticles. The hydrogen enters from distribution channels in the surrounding moderator block, 
Fig. 2.1. The control system under consideration uses a model based algorithm which relies 
on an accurate representation of the reactor in calculating the control signals. Because the 
design is unique for the application, well established thermal-hydraulic models may not be 
directly applicable. The objective of this investigation is to develop a thermal hydraulic 
model of a reactor fuel element which can, in turn, be expanded to include the entire reactor 


system and be incorporated in the controller. 


The thermal hydraulic model serves two purposes in this type of reactor control. One, 
it Is used to describe the heat transfer processes and provide the information to prevent tem- 
peratures from exceeding design limits. Two, it provides the temperature and pressure data 
required to quantify the reactivity feedback associated with the hydrogen coolant. It should 
be capable of predicting coolant flow distribution, temperature, pressure and density during 


steady state and rapid transient operations. 





1.2 METHOD 


To be an effective part of the controller, the model should function in real time in a 
stable and robust manner. These requirements place competing constraints on its form. It 
must contain sufficient detail to give an accurate representation of the fluid flow and heat 
transfer processes, but yet be simple enough that its outputs are available when required by 
the controller. Thus, constructing the model is a stagewise process. First, one must deter- 
mine what aspects of the fuel element should be included to give an accurate picture. Then, 
given the minimum requirements, the model must be tailored to meet the time constraint. 


The work described here addresses the first of these tasks. 


A detailed reference model is constructed and used to evaluate the sensitivity of the 
model predictions to various simplifications. One example is radiative heat transfer in the 
solid phase. The results of the reference case and a simpler version without the effect can be 
compared to determine if this mode of heat transfer contributes significantly to the behavior 


of the fuel element and if it must be modeled explicitly. 


The building of the reference model is the major emphasis of this work. Future devel- 
opment can then take the results and perform the simplifications required for a real time, pre- 


dictive, model. 


1.5 ISSUES 


The specific issues investigated can be summarized in two categories, hydraulic and 


heat transfer. 





Hydraulic issues include: 

a. Flow distribution: Does the model accurately reflect the flow distribution of 
hydrogen? 

b. Pressure losses: What are the major sources of pressure loss and how can they 
be modeled? 

c. Coolant flow: Can coolant flow in the particle bed be considered one dimen- 
sional or are two dimensions required? 

d. Hydrogen compressibility: Must the model explicitly represent the coolant as 
a compressible gas or are pressure changes small enough that incompressible models are ade- 


quate? 


Heat transfer issues include: 

a. Interphase heat transfer: Can heat transfer between the gas and solid phase be 
modeled as purely convective? How should the heat transfer coefficient be calculated? Is 
radiative heat transfer important? 

b. Intraphase heat transfer: Is there significant conduction and radiation between 
particles in the solid phase? How should transfer within the fuel particle be simulated? 

c. Discretization: Is a fine node analysis required or can fuel element behavior 


be represented with a few or single node model? 


The chapters that follow describe the details of the packed bed reactor and then discuss 
the development of the mathematical representations of the fuel element. These are then 
combined in a computer code to form the reference case used to evaluate the questions 


above. 
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СНАРТЕК2 
BACKGROUND AND APPROACH 


2M OBJECTIVE 


Before one can abstract a thermal-hydraulic system to a series of predictive mathemat- 
ical relations, the design and operation of the actual system must be clearly understood. The 
objective of this chapter is to present the physical details of the packed bed reactor and the 


fuel particles within it. 


Space based reactor designs have been evolving since the ROVER program of the 
1960s and its subsidiary NERVA project. After a lull in the 1970s and early 1980s, United 
States interest in the area has been rejuvenated. Major portions of the ongoing studies are 
being conducted at the Brookhaven and Sandia national laboratories. The Brookhaven work 
(L-1, H-1) analyzes and discusses reactor system designs, capabilities and design philosophy. 
Sandia's investigators have constructed prototype fuel elements and tested them in the Annu- 
lar Core Research Reactor (ACRR). The former provide the best understanding from a sys- 
tem perspective while the Sandia project gives a detailed dimensioned design suitable for 


modeling. 


2.2 THE REACTOR SYSTEM 


2.2.1 Design Philosophy 


Space based reactor systems are by necessity constrained to be lightweight and com- 


pact. This translates to a design favoring high power density, low material density and large 
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heat transfer surface per unit volume. For example, systems launched by the space shuttle 
are limited to a cylinder 4 m in diameter and 17 m long with a maximum mass of 29,500 kg. 
(H-1). Particle bed nuclear reactors are ideal for this application, because of their small size 
and their ability to deliver twice the specific impulse of a chemical rocket (H-1) in open cycle 


applications. 


Proposed systems have included both open and closed cycles (L-1). Open cycle sys- 
tems have the advantage of low weight and simplicity, since there is no need for recovery 
and reuse of the coolant. However, reactor life is then determined by the amount of coolant 
tankage. Typically these systems are considered for pulsed operation or propulsion designs. 
Closed systems have been proposed for continuous electrical power generation. Helium and 
hydrogen are the usual coolants. Hydrogen is generally preferred for the higher power 
pulsed operation systems. The pulsed open cycle variant is of principal interest in this inves- 


tigation. 


2.2.2 Pulsed Reactor Design 


Open system pulsed operation designs are discussed in references P-1, L-1 апа В-1. 
The basic cycle includes compressing liquid hydrogen, providing two stages of heating and 


directing the high temperature coolant to the power conversion unit. 


The core consists of fuel elements mounted in a moderator block which 1s in turn sur- 
rounded by a series of control drums and a reflector chamber. The entire assembly is encap- 
sulated in an aluminum pressure vessel (P-1). The core is a right cylinder .79 m in diameter 
and 1.03 m long (P-1). Other designs place the reflector outside the pressure vessel (H-1). 


The moderator block (Fig. 2.1) can be ZrH,, BeH, or LiH. Either a single drilled piece, or an 
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assembly of smaller modular blocks, the moderator is channeled to create a coolant flow path 
and sites for the fuel elements. Nineteen elements are arranged in a hexagonal array around a 
central cell (Depending on one’s perspective this may be considered a triangular arrange- 
ment). Control drums are used in lieu of control rods for cempactness and simplicity of 
operation. A portion of each drum surface is coated with a strong neutron absorber. Rotating 


the absorber away from the fuel assemblies is the equivalent of withdrawing a control rod. 


The fuel elements, Fig. 2.2, are composed of two concentrically mounted cylindrical 
retention pieces, referred to as frits. The frits are porous enough to allow coolant flow but 
strong enough to function as retainers for the fuel particles in the annular space between 
them. The space between the outer frit and the moderator block forms an inlet plenum for 
the coolant. Similarly, the region inside the inner frit forms an outlet plenum from which the 
hot exit gasses are extracted and directed to the energy conversion device, most likely a tur- 
bine, or are exhasted for propulsion. The outer cold frit is fabricated by sintering micron 
sized stainless steel particles. However, because of the high exhaust temperatures 
(approximately 2,000 K) the inner, hot frit poses a material selection problem. Drilled rhe- 


nium is used in the currently favored design. 


The fuel particles are similar to the coated fuel particles used in high temperature gas 
reactors (HTGR). Each consists of a central kernel of fuel coated with two pyrographite lay- 
ers and zirconium carbide. But unlike the HTGR the 500 micron particles are used directly 


in a packed bed configuration not combined with graphite structure to form a ball or block. 
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2.2.3 Hydraulics and Heat Transter 


Liquid hydrogen is stored in a separate container. When power operation is initiated, 
the hydrogen passes from the tank to a turbo pump where it is compressed to a supercritical 
fluid with a pressure as high as 6 MPa. It is discharged to the reactor reflector, the first heat- 
ing stage. From the reflector, the coolant passes into the reactor core. Once in the core, the 
hydrogen flows axially through coolant passages located between the fuel elements and the 
moderator block. These channels are slotted lengthwise to create a gas flow path to the inlet 
plenum throughout the entire length of the fuel element. Each coolant channel) feeds three 
fuel elements and six channels feed each fuel element as shown in Fig. 2.1. The inlet plenum 
provides communication for the coolant over the entire surface of the cold frit. Pressure 
losses in the core to this point can be considerable, since the flow has made three 90 degree 
direction changes in the core by the time it gets to the cold frit. At this point the hydrogen 
flows radially inward through the packed bed of fuel and hot frit. Collected in the outlet ple- 


num, the hot gas is exhausted to the turbine. 


Two alternatives exist for powering the turbopump (S-1,S-2 and S-3), a bleed system 
and a topping system. The bleed system takes a small portion of the reflector outlet gas, 
expands it through a turbine and exhausts it to the surroundings. The topping system 
expands all of the reflector gasses through a turbine and discharges them to the reactor core. 


The larger volume is used to reduce the pressure drop. 


One of the design imperatives of any reactor is to match coolant flow and power distri- 
bution. This matching for a packed bed reactor provides a nearly constant temperature 
throughout the length of the outlet plenum. Zoning the fuel, shaping the coolant flow, or 


combining the two methods can be used. The cold frit provides the best place to shape the 





flow distribution. Varying the porosity and/or thickness effectively controls the pressure 
drop across the frit to achieve the desired flow characteristic. The cold frit essentially serves 
as a distribution manifold. This approach is most easily applied if the flow through the 


packed bed of fuel is essentially mward (radiai not axial). 


The majority of the pressure drop in the open cycle is in the moderator block (B-1). 
Within the fuel assembly, the cold frit is the greatest resistance to flow. The hydrogen enters 
the fuel element at approximately 300 K and at 2 to 6 MPa depending on the design. At full 
power the exit gas is in excess of 2,000 K. The pressure drop through the element is on the 


order of 50 kPa. 


The packed bed reactor can operate at very high power with relatively small tempera- 
ture changes across the gas film and the fuel particle (because of the close proximity of the 
coolant and the heat deposition). This design gives a heat transfer area of 7,000 to 10,000 
m’/m’ of fuel (P-1). Powell's analysis shows that power densities of 10 GW/m’ of fuel can 
be achieved with gas film ATs of approximately 150 K (P-1). The packed bed design also 
responds well to the rapid transients associated with pulsed operation (P-1). The fuel par- 
ticles are capable of enduring repeated temperature ramps of 1,000 K/s. The reactor postu- 


lated by Sandia is expected to execute a power up ramp with a period of about 0.6 s. 


2.3 THE PIPE EXPERIMENT 


Investigators at the Sandia National Laboratory are conducting a series of experiments 
titled the Pulsed Irradiation of a Particle Bed Element (PIPE). The objective of the tests is to 
evaluate a packed bed fuel element and determine if the design is viable from the standpoints 


of temperature distributions, flow distributions and coolant/fuel compatibility (V-1). The 
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experiments place a prototype fuel element for the reactor design discussed previously in the 
Annular Core Research Reactor and pulse the reactor. This will create power densities up to 


2 GW/m! in the packed bed element. The element is full size in the radial direction; in the 


axial direction it has been reduced by two thirds. 


The PIPE fuel element is the basis for the model being developed in the work reported 


here, since it provides a complete set of consistent dimensions and performance predictions. 
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A cylinder .35 m in diameter and 4.19 m long houses the experimental apparatus, Fig. 
2.3. It consists of a self contained pressurized hydrogen circulation system. The hydrogen 
flows from blowers in the upper portion of the cylinder down around the middle section com- 
ponents to the lower section. There it enters the top of an inlet plenum fonmed by the cold 
frit and an annulus of moderator. Flow continues radially inward through the cold frit, fuel 
and hot frit before entering the central outlet plenum. This plenum discharges the hot gas to 
a heat sink contained in the middle section. The heat sink is formed by a bed containing a 
number of steel balls (originally at a low temperature). From the heat sink it passes through a 
flow meter and back to the suction side of the blowers. Parahydrogen is used, because of its 


better heat transfer properties. 


The coolant entering the inlet plenum is pressurized to 2 MPa and is at a temperature of 


300 K. Figure 2.4 shows the dimensions of the PIPE fuel element. 


The cold frit acts as a distribution manifold to match the hydrogen flow and axial 
power profiles. In this case, the frit is of uniform porosity, .685, but varies in thickness from 
1.70 to 2.36 mm to adjust the pressure drop. It is fabricated by sintering 2.5 [um diameter 
particles of 316 stainless steel. The hot frit is a uniform rhenium tube with electro-arced 


oblong holes which provide 23.3% open area for flow of the 2,000 K exit gasses. 


The fuel pellet used for the experiment is shown in Fig. 2.5. The uranium carbide fuel 
kemel is enriched to 93% and is coated with two layers of pyrographite. The inner one 1s 
low density to accommodate fission products. An outer coating of zirconium carbide, which 


was developed to improve fuel coolant compatibility in the NERVA project, is also present. 
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2.4 EXISTING MODELS 


Several models have been developed for annular fuel elements. One early model was 
formulated in Britain for the design of a gas-cooled fast reactor (A-2, B-2). More recently, 
Brookhaven researchers (B-1, L-1)(in conjunction with their designs for space reactors) have 


modeled packed bed elements similar to those being tested at Sandia. 


Each of these models was targeted to a specific purpose which in tum influenced the 
assumptions made and the detail included. The British model was concerned only with flow 
distribution calculations and did not address the heat transfer issues. One of the Brookhaven 
models coupled the neutronics and heat transfer to create the control codes KINETIC and 
SPHEAT (L-1). In these, the fuel element is lee radially and axially to allow for vari- 
ation in the power distribution. However, hydraulics was handled by assuming that the cool- 
ant flow matched the power in each axial slice. ‘The thermal-hydraulic code described in 
references (A-1) and (B-1) uses the same partitioning scheme, effectively giving a 2 
dimensional model, but again assumes the flow in the fuel matches the power distribution 


and uses this to back calculate the required pressure drop characteristic of the cold frit. 


The model developed in this investigation is dedicated only to the thermal hydraulics. 
Its interaction with the neutronic model is limited to passing results such as hydrogen temper- 
atures, pressures and densities. With this in mind, the intention 1s to build a more general 
model that operates using only input information that can be supplied by system 
instrumentation (the inlet coolant temperature, the inlet pressure and the discharge pressure). 
In this way, one can determine if the flow and power match for actual operating conditions 
and not just for expected design conditions. The model is also general enough to allow axial 


as well as radial flow and to include heat transfer between the particles of the solid phase. 
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25 MODEL STRUCTURE 


2.5.1 The Premise 


The naturai tendency in modeling ts break the problem down, look at the individual 
pieces and then reassemble them. In the case of the packed bed fuel element, two subdivi- 
sion methods are available, separation by component or by phase (solid or gas). These can 
be combined in several ways. The most detailed is to consider each phase in each 
component. But this makes it difficult to put the components back together. Another option 
works with each phase in its entirety and then joins them. The latter lends itself to the appli- 


cation of the fundamental laws of physics and easy model unification. 


The model developed in this investigation treats each phase independently and joins 
them at the common point, solid to gas heat transfer. Thus one could look at the fuel element 
as a large heat exchanger. The gas phase represents the shell side and the solid phase with its 
heat deposition is the tube side. Heat flow in the gas phase can be convective and/or conduc- 
tive and is independent of the solid phase except at the connection point. Similarly, heat flow 
by conduction and radiation in the solid phase is only influenced by the gas phase at the point 
of heat transfer. This representation is one of the basic tenets of the Momentum Integral Net- 


work method developed by Van Tuyle (V-2, V-3). 


The model is based on the premise that the fuel element must obey the fundamental 
principles described by mass, momentum and energy balances. To facilitate the application 
of the balances the element is divided into an array of control volumes. The volumes are 
arranged to coincide with the components of the fuel element. The physical features of the 


components determine the way the terms of the balance equations are computed in their 
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respective volumes. Thus, differences between model components can be easily incorpo- 
rated. For example in the momentum balance, a Darcy friction factor for pipe flow is used to 
calculate the frictional losses in the plenum control volumes, and an Ergun relation for 
packed beds is used in the volumes containing fuel particles. The matrix of control volumes 


and the continuity of conditions at the interfaces provide the required unifying structure. 


The gas phase model uses all of the balances: mass, radial momentum, axial momen- 
tum, and energy. The inclusion of these balances accounts for the coolant compressibility. 
However, since the solid phase is fixed and virtually incompressible, only a heat balance is 
required. For transient analysis, each phase is advanced through the time steps in parallel. 
The connection is made in the source (sink) term of the energy balances. Using the new time 
step solid and gas temperatures, the heat transfer is calculated implicitly and used as the 


source (sink) term for the current time step calculation. 


Although the basic equations are easily derived and manipulated in continuous form, an 
analytical solution for the fuel element would be very difficult. Therefore, the balances are 
discretized in a finite difference form based on the control volumes. This is readily adaptable 


to computer solution. 


2.5.2 The Staggered Grid 


In flow problems, the discretization process is greatly enhanced by defining the vari- 
ables at different points. For example, velocities are defined at points between points at 
which pressures are defined. The resulting arrangement is often referred to as a staggered 
grid. Its advantage lies in the fact that the variables are defined where they are needed. In 


the momentum balances, the pressure difference across the velocity control volume provides 
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the source term. Using the staggered grid, pressure is defined on all faces of the velocity 
control volume. An aligned grid would require an estimation of the face pressures. This fre- 
quently introduces large errors in the numerical solution (P-2). The same is true for the mass 
balance where the control volumes are centered on the pressure nodes. The balance requires 
the mass fluxes crossing the volume faces. The staggered grid provides a velocity at each 


face for computing the required mass flux. 


Figure 2.6 shows the staggered grid used in this model and the control volumes used in 


each of the fundamental balances. Cylindrical coordinates are be used throughout the model. 
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2.6 SUMMARY 


Proposed space based reactor designs are centered on the packed bed fuel element with 
its high power density. Although the entire reactor has never been built, Sandia National 
Laboratory has constructed and tested annular particle fuel elements. These prototypes are 


the basis for the modeling and the analysis performed during this investigation. 


The model determines the temperature, pressure, density and velocity distributions for 


both the gas and solid phases through the solution of discretized mass, momentum and 


energy balances applied to the staggered grid of control volumes shown in Figure 2.6. 
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CHARLER3 
THE GAS PHASE 


3.1 OBJECTIVE 


The gas phase is the most complicated part of the model, since it contains both hydrau- 


lic and heat transfer considerations. This chapter describes the underlying assumptions and 


modeling of the gas phase. The discussion uses the continuous form of the equations. The 


discretized form and solution method will be presented in chapter 4. 


3.2 CHARACTERISTICS INCLUDED IN THE MODEL 


The model of the packed bed fuel element should include all of the features that signifi- 


cantly affect the hydraulic or heat transfer response of the system. Additionally, it should be 


able to track transients as well as solve for the steady state condition. So, the tume dependent 


terms must be considered. The primary physical features of interest are: 


d. 


b. 


f. 


8: 


Friction effects 


Manifold “action” of the plenums 


. True two dimensional flow 
. Pressure losses due to changes in direction 


. Spatial acceleration due to reduced flow area and due to density changes 


Gas film temperature differences 


Enthalpy changes due to pressure changes and heat transfer 


In order to incorporate these aspects in the model, the basic equation for each of the balances, 


momentum, mass and energy is developed. Then each term is reviewed and modeled includ- 


ing these characteristics when appropriate. The description of the model details that follow 
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uses a similar approach. However, items strongly affected by the geometry and 
approximations used in the solution method are not specifically addressed until the model is 


discretized. Spatial acceleration is one example. 


Several assumptions are made to keep the problem tractable and eliminate insignificant 
effects. Given the fuel elements inherent symmetry, it is assumed that the entire problem can 


be treated as cylindrically symmetric. Thus, only the radial and axial variations are included. 


Gas velocities are expected to range from one to several hundred meters per second. 
As a result, the primary heat transfer mechanism in the gas phase is convection. Conduction 
and radiative processes are only considered significant in the solid phase. The model con- 
tains other assumptions that are specific to certain aspects. These are discussed when 


encountered in the model development. 


Boundary conditions are also factored into the modeling process. The model assumes 
that only the inlet pressure, inlet temperature and outlet pressure are known in addition to the 


physical dimensions and the initial temperature, pressure and velocity distributions. 


3.3 THE MOMENTUM BALANCES 


3.3.1 The General Form of the Equation 


The model contains two separate momentum balances, one for radial flow and one for 
axial flow. The development and application of the two is identical; the only difference 
being the direction of the velocity field. The momentum balance discussed below is the gen- 


eral one which can be applied to either the radial or axial velocity field. 
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The general form of the balance equation is based on a lumped parameter fixed control 


volume approach. It is derived from the General Transport Equation, Eq. 3.1 (T-1). 


D — d — — --» —+ 
= f ) e 
od | Jeena | | jen +} (r,t)v, * ndS Sl 


f(r,t) is a general transport function and v, is the velocity of the fluid relative to the control 


volume. In this case, since the volume is fixed, it is the just the velocity of the fluid. V.S and 
t represent the volume, surface area and time respectively. This equation can then be made 
specific to the momentum balance by substituting pv , the linear momentum per unit volume, 


for the transport function and realizing that 
я =- | pv, * nds 3.2 
S 


where m 1$ the net mass flow rate into the control volume and » is the mass in the control 


volume. The result is: 


Duce: ud. 5 e 
DUE (En д - 2 ту, 9:3 


where ii; is inward the mass flux through the ith face of the control volume (cv) and v; is its 


velocity component in the direction of the momentum balance. The product then represents 


the momentum flux through the i^ face. According to Newton's second law the time rate of 
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change in momentum is equal to the sum of the applied forces. Applying this to Eq. 3.3 
gives the final form of the momentum equation (This 1s a synopsis of the more detailed deri- 


vation given in reference T-1). 


In other words, the rate of change of momentum within the control volume 15 equal to the 
sum of the forces acting on the control volume and the net change in the contained momen- 
tum due to flow through the boundaries. Therefore, the effects to be included in the model 


fall into 2 general categories the source terms (pressure, friction) and the flux term. 


The control volume used in the momentum calculations is centered on the appropriate 


velocity as shown ın Fig. 2.6. 


3.3.2 The Source Term 


The source term accounts for two effects, friction and pressure. Because of the geome- 
try and the variation in the structure of the fuel element components, the evaluation of this 


term is dependent on the control volume and the velocity variable involved. 
The most easily analyzed force in the source term is the effect of pressure. The net 


force in a given direction is simply the pressure difference between opposing faces multiplied 


by the effective area. The staggered grid arrangement defines pressure on each of the faces, 
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so estimation is not required. However, geometry creates some pitfalls that must be under- 
stood to prevent introducing errors. Pressure in particular is difficult to work with in cylin- 
drical coordinates and still maintain the proper vector relationships. Therefore, rather than 
resolve the pressure forces to an x-y coordinate system, the fuel element is treated as though 
it were slit axially and then flattened like a sheet. This way all radial forces act in the same 
direction. Care must be taken to ensure that this treatment does not distort other physical 


features such as control volume face areas. 


The control volumes of the gas and the solid phase are intimately intertwined. The 
dimensions of the volume include the space occupied by the gas and the solid. The volume 
of the gas and its associated control volume are related by the void fraction, € , a geometric 


parameter of the system. 


КЕ bee, ae 3.5b 


The same relationship also holds for the area of a control volume face and the effective area 
of the gas. It must be remembered that in a momentum balance on the gas phase that the 


pressure difference of interest acts on the gas area not the combined gas solid area. 


This relationship is also important when velocities are being considered. This model 
works with the so called "superficial velocity" which is the velocity the gas would have if the 
control volume were completely empty. Again, superficial velocity 1s just the product of the 
void fraction and the velocity of the gas in the spaces between the solid phase, the interstitial 


velocity. 
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The geometrical arrangement of the fuel cell and coordinate system also means that the 
differential control volumes shown in Fig. 3.1a are not rectangular parallelepipeds. In a zero 
gradient pressure field, a net force would appear to be exerted on the volume if only the top 
and bottom faces were included in the area calculation. This is not physically correct as a 
zero pressure gradient should yield zero for the pressure force. The model remedies this by 


including the side faces in the following manner. 


Given the pressure diagram, Fig 3.1b, for a two dimensional differential volume, the 
area of the top and bottom faces is approximated by the chords T and B respectively. The 


downward direction force balance is: 
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Figure 3.1a A Typical Control Volume 





Figure 3.1b Control Volume and Pressure Field 
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bottom 
This can be reduced to 


Net Force=P M+P 


top bottom 


M 7 


where M is the mean chord length, (T +8)/2. This is not a problem for axial pressure forces, 


since the opposing faces are of equal area and the side faces are parallel to the direction of 


pressure application. 


Two mathematical representations of the friction forces are used, depending on the 
location of the control volume being considered. The flow in the inlet and outlet plenums is 
viewed as pipe flow with friction forces proportional to a friction factor. ‘The frits and the 
fuel particles are treated as packed beds using the pressure gradient described by the Ergun 


relation. 


The inlet and outlet plenum pressure drops due to friction are given by the relation 


L \pv’ 
sr š 


where: 


P = Pressure 


ч, 
І! 


Friction factor (Darcy) 
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p - Density 
v = Velocity 


= Length over which the pressure drop is experienced 


D = Equivalent pipe diameter 


Once the pressure drop is known, the frictional force can be determined by multiplying by 


the appropriate area. 


All of the values required to calculate the friction force are either known from the 
geometry of the control volume and the previous time step values of the unknowns or are 
variables in the final solution. The exception is the friction factor which requires an interme- 


diate computation. 


Blasius and McAdams proposed two of the more common correlations for the friction 


factor of smooth pipe. 


f2-0484Re ^? (McAdams) 3.9a 


f2-0316Re ^? (Blasius) 3.9b 


The McAdams is generally used at Reynolds number, Re, greater than 30,000 and the Blasius 
at lower values of Re. However, direct use of these relations can be improved upon. The 
Moody diagram shows friction factors as a function of the pipe wall roughness as well as the 


Re. As pipe wall roughness increases the exponent becomes less negative. 
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In an effort to preserve the easily computed form of the McAdams and Blasius equa- 
tions and account for the roughness of the frits, a new correlation line was fit to the log-iog 


based Moody plot. The result was. 


f=0.138Re” 3.10 


This revised correlation also has the advantage that it is applicable over a wider range of Re 
values (10° to 10°), thus preventing the potential problems of discontinuities in the friction 
term associated with the standard application of the Blasius and McAdams formulas. For Re 
less than 1,000, a constant value of 0.0482 was assumed for the friction factor. A plenum Re 


this low should only be encountered near zero flow initial conditions. 


The Colebrook equation (C-1) would also have worked and does include the surface 
roughness dependence. But its transcendental form makes it difficult to compute. Although 
the friction model for the plenums appears crude, it is probably as accurate as required. The 
numerical method will use previous time step velocities to calculate the Re as friction losses 


in the plenums are small compared to those in the frits and packed bed. 


There are two ways to look at flow in packed beds, either as flow through tubes with 
irregular cross sections or as flow around objects. The former approach has generally been 
more successful (B-3). The Ergun relation, Eq. 3.11, is a widely used example of the tube 


model (E-1). 


АР 150pvo(1 Y, 1.75pvol vi (1 - £) T 
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where: 


AP/L = Pressure drop per unit length 
D, - Particle diameter 
vo = Component of the superficial velocity in direction of the momentum 
balance 
|v| = Superficial velocity magnitude 


This relation is really the smooth blending of two correlations, the Blake-Kozeny for low Re 
and the Burke-Plummer for high Re. So in essence, the first term is a laminar term and the 
second is a turbulent term (B-3). In extremes of flow, one of the terms will dominate. As 
with the plenum calculations, the force on the control volume is determined by multiplying 


the pressure difference by the appropriate area. 


Several assumptions are implicit with the use of Eq. 3.1]. First the cross sectional area 
for flow is assumed to be constant. This is approximated in the model by using the mean 
cross sectional area of the control volume for radial flow. The condition is satisfied exactly 
for axial flow. Second, the packing is uniform everywhere in the fuel element. This removes 


the effects of channeling (B-3). 


Third, the ratio of the particle diameter to the effective bed diameter is small (B-3). 
Thus wall bypass effects are assumed to be negligible. These effects arise from the variation 
in void fraction as a function of radial position. It is a maximum of | at the wall where there 
is only a single point of contact between the wall and each sphere. The minimum of 0.2 
occurs half a particle diameter in from the wall. The void fraction then shows a damped 


oscillatory variation as one moves radially inward and eventually reaches the mean value. 
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The unity void fraction at the wall essentially creates a low resistance channel which allows 
flow to bypass the bed. Reducing the particle diameter to bed diameter ratio minimizes this 


effect. 


The Ergun equation uses the superficial velocity. In order to preserve the vector nature 
of the calculation, the velocity component in the direction of the balance is used. When 
velocity is squared in the second term, the magnitude of the velocity vector as well as the 


velocity component along the direction of the balance are used. 


The Ergun relation forms the basis of the model’s friction calculations in the frits and 
fuel particle bed. All of the terms in the Eq. 3.11 are known either from the geometry or the 
previous time step. The method of solution discussion will address how the gradient is com- 


puted given this basic relation. 


Achenbach (A-2) also proposed a relation to predict packed bed pressure drops which 
was adopted by the Brookhaven researchers. It gave 25% higher gradients than the Ergun 
relation and was chosen as the conservative method (B-1). This model uses the Ergun which 


has worked well in the Sandia PIPE experiment preparations. 


3.3.3 The Flux Term 


The net momentum flux into a control volume is determined by looking at the mass 
flow rate through each face and its velocity then combining them vectorialy. The mass flow 


rate is computed by equation 3.12. 


4] 





т, =руа ag 


mv.-ovglv; 5.126 


In this case, the area, a, 1s the overall area of the control volume face. Note also, that v, is 


the component of the superficial velocity perpendicular to the face of the control volume. ri 
is positive for flow into the volume and negative for outward flow. Since velocity is a vector 


with a sign convention, care must be taken when determining the sign for this term. 


Once the mass flow rate is known, multiplying by the velocity component in the direc- 
tion of the momentum balance, v;, yields the momentum flux for the face, Eq 3.12b'. Figure 
3.2 illustrates the momentum flux for a typical control volume and a radial momentum 


balance. 


Figure 3.2 shows the use of the staggered grid and its complications. Fluids enter the 
West face with two different velocities and densities depending on which control volume 
they enter from. No attempt is made to model a more detailed gradation of density and 
velocity, since the lumped parameter approach assumes the value of a variable to be the same 
everywhere in the volume. All variables in Eq. 3.12b are known or are a variable in the final 


set of solution equations. 


| v. as used in the model development and results that follow is defined to be the superfi- 
cial velocity. This is incorrect. The formulation of the momentum balance used in the model 
requires that v; be the interstitial velocity. As a result, the model evaluates the time 
derivative and momemtum flux terms incorrectly. Refer to section 7.5, Recommendations. 
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3.3.4 The Manifold Effect 


One of the principal functions of the inlet plenum and cold frit is to serve as a flow 
distribution manifold. Similarly, the outlet plenum serves as a collection manifold. This 
raises the question ot whether or not the general momentum balance accurately predicts the 
radial flow as a function of position and whether or not improvements can be made. This is 


especially important given the requirement to match the power and flow distributions. 


Rephrasing the question for the inlet plenum, is the average axial component of 
momentum of the radial (lateral) stream leaving the control volume in Fig 3.3 the same as 
that of the fluid at the center of the volume where the predominant flow is axial? The model 
as constructed to this point assumes that the averages are the same. In other words, the axial 
velocity component of exiting fluid, v,, is equal to v, in Fig 3.3. The other extreme says v, is 
O and the exiting stream carries no axial momentum. Bajura noted this issue and included it 


in his models for piping distribution manifolds (B-4,5 and 6). 


Bajura derived an analytical expression for manifold flow based on an axial momentum 
balance which included factors to account for deviations from ideal behavior. He then per- 
formed a series of experiments to evaluate the constants. He assumed steady state conditions 


and started with the following balance for the dividing header volume shown in Fig 3.3: 


44 





Axial Stream 


C 
Е 








= a a a a a m - a œ - a- ea aa a aa n+ 


Radia! Stream 


Figure 3.3 
Manifold Control Volume 


(Adapted from 8-4, B-5, B-8) 
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| P,dA -$ P,dA +} т„ЧА =} pv;dA -$ ру, ЧА + | pv,v,dA Du 
Ai А, AZ A, A, va 


The Left hand side contains the source terms and the right hand side the flux terms. 
Then he defined two correction factors that involve an area weighting and the relationship 


between v, and v,. 





ajo 24А 3.14a 
= 


l 
=} vv dA 3.145 
3 А, i 


моу з: 


B and y are the axial momentum correction factors for the axial and radial streams respec- 


tively. B Accounts for any deviations from a flat flux profile in the axial stream and y is a 
measure of the axial momentum of the radial flow stream as a function of the axial 
momentum per unit volume in the control volume. Substituting these relations, applying a 
continuity constraint and assuming that conditions on opposite sides of the control volume 
are related by a first order Taylor series expansion, Bajura’s balance becomes after consider- 


able algebra: 


ЗоигсеТепи$ = -ру А, бх E — (2B — y)vopA, à 3.15 


d : . i 
T accounts for changes in the flux profile due to plenum entrance effects as a function 


of position. D becomes constant within the first 20% of the header length when fully devel- 
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oped flow is established. It is very dependent on the header geometry. (2B — y) 1s the overall 
momentum correction factor, 0. If v, equals vo and the velocity profile is flat across Aj, then 
Y equals 1. If the velocity profile across A, is flat, then D equals 1. A similar analysis can 
also be made for the outlet header. This particular condition is exactly that contained in the 
model if the momentum correction factors were not included. Bajura’s experiments found 
that for fully developed incompressible flow the values of Y were 0.95 and -0.66 for the inlet 
and outlet headers respectively (B-4). The corresponding values of 0 were 1.05 and 2.66 


(B-6). These are consistent with plug flow and a D of 1. 


The fact that y for a dividing header is approximately .95 says that v, is less than v; 


(B-4). This can be explained if one assumes that the lower kinetic energy boundary layer is 
preferentially turned to radial flow by the virtue of its position relative to the bulk axial flow 
and the radial flow openings. y for a combining, outlet, header tends to -.66 and represents a 
fundamentally different interaction between the radial and the axial streams. Ín a combining 
header, the radial stream carries virtually no axial momentum into the plenum control vol- 
ume. But it does add to the mass in the header and effectively reduces the flow area avail- 
able to axial flow stream entering the volume at the upstream face. This necessitates a 
velocity increase and a further pressure drop in addition to that required to accelerate the 
ЕЕЕ flow and counter friction effects, thus the large value of 9 in a combining header 
(B-4). y should not be viewed as implying the axial momentum of the radial flow stream 15 


in the opposite direction of the plenum flow and momentum. 


Datta and Majumdar (D-1, M-1) took Bajura's results one step further. They trans- 


formed the equation to a finite difference form and relaxed the requirements for constant 
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porosity and lateral (radial) flow resistance. Using this and Bajura’s recommended values for 
d . ` : 
8 they found that the entrance effect term, = did not significantly affect the results and 


assumed [ was constant for a given header configuration and flow. 


The modei of the fuel element combines the preceding formulations and adds the time 
dependency to transform the general momentum equation, 3.4, to the following form for the 


inlet and outlet plenums. 


^ NR = > 
FIM ze Enc af pv,dA -$ pv dA -$ oaa 3416 
dr са К deu А, А, Ay 


Refer to Fig 3.3 for notation. The radial momentum balance for the plenums retains the form 


of Eq. 3.4. 


Although this model is based on incompressible flow, it still represents an improve- 
ment to the unmodified momentum balance. Anticipated gas velocities do not exceed 30 to 
40% of the speed of sound, so the compressive effects should be minimal. The value of B is 
assumed to be 1. However, if solutions indicate flows are near the laminar turbulent trans- 


ition B can be adjusted to achieve a better fit of predictions to the data. 


From the designer’s perspective, the two variables with the greatest influence on flow 
distribution are the ratio of axial to total radial flow area and the resistance to radial flow. 
Ideally, the larger the radial flow resistance and the smaller the area ratio the easier it is to 


achieve the desired distribution of flow (B-6). The fuel element uses the radial flow 
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resistance of the cold frit to control the flow distribution. The cold frit pressure drop is over 
10 times the pressure drop of the rest of the fuel element. This high pressure drop is the price 


of maintaining flow distribution stability under a wide range of operating conditions. 
3.4 THE MASS BALANCE 


3.4.1 The General Form of the Equation 


The mass balance accounts for the amount of hydrogen in the system and ensures that 


in steady state the mass entering the system equals the mass leaving of the system. 
The general equation describing the mass conservation 1s derived from the General 


Transport Equation, Eq. 3.1. In this case the transport function, f(r, t) is replaced by p, the 


density. The equation becomes: 


Dm d | f ae 
—— === dV + “nds 317 
реа и!) 2 Hi j 


Substituting Eq. 3.3 this reduces to: 


I 
br (| = > m, 3.18 
Dr a 


Since the substantial derivative of the mass equals zero, the continuity balance can be stated 


as: 
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The time rate of change of the mass in the control volume is equal to the net mass flux into 


the volume. 


There are no source terms in the equation. This is consistent with the absence of any 


chemical or nuclear reactions that produce a significant amount of hydrogen. 


The continuity balance could theoretically be used on any arbitrarily defined control 
volume. However, for the convenience of the solution algorithm, it will be applied to a pres- 


sure centered control volume as shown in Fig 2.6. 


3.4.2 The Flux Term 


As done with the momentum balances, the mass flux across a control volume surface 
is the product of the density, velocity perpendicular to the face and the face area, Eq. 3.12. 
Since it is the superficial velocity, the total face area should be used. The difference from the 
momentum case is that the velocities are defined on each face of the control volume in accor- 
dance with the staggered grid simplifying the evaluation of the flux terms. Thus, everything 
required is known from geometry or is a defined variable obtained when the overall solution 


is completed. 


50 





3.5 THE ENERGY BALANCE 


3.5.1 The General Form of the Equation 


The energy balance serves two purposes. First it provides the additional equation 
required to solve for the state variables of the hydrogen. The momentum and mass balances 
are sufficient to obtain a hydraulic solution for an adiabatic incompressible flow system. 
However, for compressible flows, since density is a function of temperature and pressure, the 
third equation in conjunction with the equation of state is necessary to obtain velocity, tem- 
perature, density and pressure. Second the enthalpy balance ts the primary heat transfer rela- 
tionship representing the connection between the solid, heat producing, phase and the coolant 


gas phase. 
The energy balance is also derived from the General Transport Equation, Eq. 3.1. This 


time the transport function is ph, where h is the specific enthalpy. Substituting and integrat- 


ing Eq. 3.1 yields 


І 
ош = 2 — > mh 3.20 
Dt dt cu { та} 


The substantial derivative (follows motion of a given mass) of the total enthalpy, H, is: 


DH dP 
---+ О +V 
Dr dt 
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where Q, the heat input, and V —, the change in the pressure volume product, are the source 


terms. Substituting in Eq. 3.20 and rearranging gives: 


d P i 
—H = +V — + n] Е 
4 | О 5 2 mh 2:22 
Thus, the time rate of change of the enthalpy in the a fixed control volume is the sum of the 


heat input, change in the PV product and the net enthalpy flux into the volume. 


The staggered grid defines all the state variables at the same points. Therefore, 
enthalpy, pressure, temperature and density are co-located and the control volume is identical 


to that used for the continuity equation. 


3.5.2 The Source Term 


There are two "sources" of enthalpy within the control volume. One is the heat gener- 
ated in the solid and transferred to gas phase. The other is the change in control volume pres- 


sure. 
The pressure dependence is the simpler to analyze. Since the control volume size is 


constant, the change in the enthalpy is the product of the volume and the change in pressure. 


The change is the pressure increment since the last time the enthalpy was evaluated. 
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When the fuel element is in operation, the larger source will be the energy transfer from 
the solid phase. As with most transport phenomena, the heat transferred to the hydrogen is 
proportional to a driving force. In this case, the driving force is the temperature difference 


between the fuel surface and the bulk coolant temperature. 


О = hA,V(T, ES Тс) 3.23 


where Q is the total energy transferred per unit time in the control volume. h is the heat 


transfer coefficient and Ay is the total particle surface area per unit volume of the bed. 


The heat transfer coefficient is a function of the mass flow rate, the shape of the par- 
ticles and the physical properties of the hydrogen. Empirical correlations have proven most 


effective for determining h. The relation used in this model is (B-3) 





Jy -0.91Re ^ y(Re « 50) 3.24a 
ін -0.61Ке ^" w(Re » 50) 3.24b 
С 2/3 
Ju A 0] 3.24с 
С, 
Ве = 3.244 
анду 


ің = Colburn j factor 
on = Heat capacity at constant pressure of the gas at bulk temperature 
: of = Heat capacity at constant pressure of the gas at film temperature 


25 





К = Thermal conductivity of the gas film 


п, = Viscosity of the gas film 
Go = Superficial mass velocity 
w = Particle shape factor, 1.0 for spheres 


The film is the boundary layer of coolant next to the surface of the fuel particle. The 
correlation approximates the film temperature by averaging the solid surface temperature and 
the bulk gas temperature. The Reynolds number used in this correlation is different than the 
one used in the momentum equation pipe flow friction term. Here, the length term is the 
inverse of the particle area per unit volume as opposed to the equivalent diameter for flow. 


The result 15 that the Re should range from 5 to 20 in the packed bed. 


3.5.3 The Flux Term 


The flux through each face of a control volume is the product of the density, superficial 
velocity, area and the specific enthalpy. These are added vectorialy to determine the net 


enthalpy flux. The final form of the balance is: 


d аР. 
EH =hAJV(Ts-T;)+V-— + 2 hh, 9:25 


i=l 


All of the terms in the balance can now be either evaluated or expressed as the product 
of a constant coefficient and a variable that will be determined when the balances are solved. 
The physical properties are generally evaluated using previous time step values for the state 


variables. 
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3.6 SUMMARY 


The gas phase model of the packed bed fuel element consists of four lumped parameter 
conservation equations simultaneously applied to an array of control volumes. The balance 


equations are shown in Table 3.1. 


Table 3.1 


The Gas Phase Conservation Equations 


= I 
mic cg 2 ту 


axial,t 


I 
+ mv 
¡al 


friction radial,t 


Note 1: 8: Inlet Plenum 0.95, Outlet Plenum 2.66, Otherwise 1.0 





The basic equations are the same for all of the volumes. The differences between com- 
ponents are reflected in the way the individual source terms are evaluated and the dimensions 


of the control volumes. 
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CHAPTER 4 
SOLUTION METHOD FOR THE GAS PHASE 
THERMAL HYDRAULICS 


4.1 OBJECTIVE 


The fundamental gas phase balances described in chapter 3 contain sufficient informa- 
tion to solve for the hydrogen temperature, pressure, velocity and density as a function of 
position when they are combined with appropriate boundary conditions and fuel element 
geometry. Many numerical schemes have been developed for the solution of fluid dynamics 
problems. This chapter describes the method chosen for this analysis, including the required 
transformation of the equations to finite difference form and the approximations imposed on 


the model. 


4.2 SOLUTION METHODS 
4.2.] Choice of a Numerical Method 


The basic problem in computational fluid dynamics is to reduce a highly nonlinear set 
of equations to a set of simultaneous linear equations that still retain the properties of the 
originals and can be solved using standard matrix techniques. All methods which are applied 
to digital computers require transformation of the continuous equations to a discrete form, 
finite difference in this model. In the discrete form, the variables can not vary continuously. 
Rather the values are defined at specific points which are paired with a control volume. Con- 
sistent with the lumped parameter approach, the value of the variable is assumed to be the 


same everywhere in the associated control volume. 
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The discretization is also the step that performs the linearization and determines how 
tme is handled. Implicit differencing uses the unknown values of the next time step and the 
current time step values; where, explicit differencing uses only current and past values to 
advance to the next time step. Implicit equations require greater computational effort than 


explicit ones. 


Solution of the discretized equations can be accomplished by iterative, predictor- 
corrector or shooting techniques. Each type has its strengths. In the end, the choice of a 


method is a trade off between stability, accuracy and computing time. 


Three numerical methods commonly found in the solution of heat transfer and fluid 
flow problems are the Semi-Implicit Method for Pressure Linked Equations (SIMPLE)(P- 
2,P-3 and P-4). The Pressure Implicit with Splitting of Operators method (PISO) (I-1,I-2) 
and the Implicit Continuous-fluid Eulerian method (ICE) (H-3). All of these methods use the 
same basic balance equations but manipulate them differently to achieve the final solution. 


Some degree of implicitness is included in all of them for numerical stability. 


The SIMPLE algorithm is an iterative technique developed by S. V. Patankar and his 
co-workers (P-2,P-3,P-4). In this method, one first guesses the pressure field values at the 
advance time, n+1. Then, the corresponding velocities are calculated using the momentum 
equations. The mass and energy balances are used to correct the pressure and velocity fields. 
The new pressure is compared to the initial guess and iteration continues until the desired 


degree of convergence is obtained. 
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The PISO method is similar to SIMPLE except that it starts with the current pressure 
field and is noniterative. The technique uses the mass balance and a predictor corrector for- 


mat to remove the linkage between the pressure and velocity equations. 


ICE is also noniterative and has been employed successfully for the solution of two 


phase flow problems in the THERMIT code (K-1). 


The PISO algorithm developed by R. I. Issa (I-1,I-2) has been adopted for this model of 
a packed bed fuel element. Since it is noniterative, the solution time is expected to be the 
shortest giving the greatest potential for use in a real time digital controller. The semi- 
implicit structure of the finite difference equations provides the necessary stability. The ICE 
method would also be acceptable. Although PISO does not allow the arbitrary accuracy 
available in iterative procedures, it is considered as accurate the differencing scheme used to 


discretize the equations (I-1). 


Unlike SIMPLE and related methods which use density corrections or update density at 
the end of an iteration, PISO forces the final pressure and velocity fields to satisfy the 
momentum and mass balances simultaneously. This feature is the key to avoiding the need 
for iteration (I-1). For a compressible flow test problem, Issa found that PISO required 19% 
of the computing time of SIMPLE and that the method was stable over a much wider range 


of time step size (I-2). 
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4.2.2 The PISO Algorithm 


The following summary of the PISO procedure is presented to motivate the discretiza- 
tion of the balance equations. Sections 4.6 and 4.7 contain the detailed sequencing of the 


equations and the boundary conditions. 


PISO uses two predictor corrector stages to advance the values of the variables from 
time n to n+1 through a series of ever more refined approximations. The first predictor uses 
the momentum balances to make the first estimate of the n+1 velocity field. It uses the cur- 


rent (already computed) tume n values for all variables except the desired velocity. 


The radial and axial velocities calculated in the first predictor satisfy the momentum 
balance but usually not the mass balance. The first corrector applies the continuity equation 
to the momentum balances. The resulting pressure equation is used to provide the first 
approximation to the new pressure and density. The second estimate of the velocities ıs also 


obtained. 


The second predictor is accomplished using the energy conservation equation. In this 
case, the enthalpy balance is used. It advances the temperature and serves as the heat transfer 
link to the solid phase. The new temperature allows one to update the equation of state rela- 


tionship between temperature, density and pressure. 


The second corrector is similar to the first corrector. It is derived by combining the 
momentum equations, mass balance and the new equation of state. The result is the second 
determination of the pressure and density at n+1 and the third approximation of the velocity 


field. 
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If one desired, a third predictor corrector stage could be added. The method discussed 
here uses only two stages. Thus, to advance the variables one time step the method makes 
three successive estimates of the velocities, two updates of the pressure and density and one 
approximation for the new temperature field. To see how this scheme works one must first 


look at the discretized form of the equations. 


4.3 DISCRETIZATION OF THE MOMENTUM AND ENERGY BAL- 
ANCES 


4.3.1 The General Equation 


The built in stability of the PISO method starts with the transformation of the balance 
equations to an implicit finite difference form which is amenable to a computer based numer- 
ical solution. The equations are arranged to advance the values of the variabies from the cur- 


rent time step, n, to the n+1 time step. 


A review of Table 3.1 shows a useful similarity in the momentum and enthalpy conser- 
vation equations. All of them equate the time rate of change of a quantity in a control vol- 
ume to the sum of the sources and the fluxes crossing the control volume surfaces. 
Therefore, discretizing a general equation using the dummy variable 6$ essentially discretizes 
all of the equations. The radial and axial velocities or the specific enthalpy can be substituted 


for 6. Thus, one starts the transformation with the general conservation equation 


d 1 
E » = source + > mb 4.1 
cv ier 
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The time derivative can be represented in discrete form by the following implicit finite 


difference 


d оф риф" E a" t! аф” 
r 90, = тл с шл е sy 42 
2 | | Št | P | 5 


The time n density is used, since the density at n+1 is not determined in the first predictor 


step of the PISO method. 


The flux terms are the product of the mass flow rate across the control surface and the 


quantity 6. The mass flow is evaluated using the expression 
mass flow rate = pva 4.3 


where у 15 the velocity normal to the control volume surface of area a . For the control vol- 


ume shown in Fig. 4.1 the $ flow rates across the control volume surfaces are 


(P v. aa), 4.4a 
(p v. aa 9) 4.4b 
(p^ v, d 0), 4.4c 
(p^ v, uua È), 444 


where the lower case letters identify the surfaces. 
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Figure 4.1 


The Control Volume for 


the General Variable Phi 
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A complication occurs in the momentum equations when $ represents a velocity. This 


results in making the flux term a function of the square of the velocity and hence nonlinear. 
Fully implicit differencing would have all of the Os, except the one in the time derivative, 
evaluated at time n+1. In order to obtain a linear equation a compromise is made. The 
velocity in the mass flow calculation is evaluated at the current time, n, and ó is evaluated at 
the n*1 time step. This also resolves another problem. If the axial momentum balance is 
performed before the radial balance, v.j,., is required to compute the flow through the north 
and south faces of the control volume but is not known. The compromise allows the use of 
the current time value of the radial velocity. The same situation exists in the radial balance. 
As a result, the flux term remains a function of velocity squared yet linear in terms of +" ^". 


This approximation is made in the flow calculations in both the momentum and enthalpy 


equations. 


The question now becomes how to evaluate a variable if it 15 not defined at the control 
volume interface, for example 6 in Fig 4.1. One scheme would be to average the values on 
either side the interface. This works in some instances and is done to estimate the current 
velocity in the mass flow term. However, this often breaks down when convection due to 
fluid flow is present. If the gas phase flows from W to O, the value of the gas density at the 
interface much more nearly resembles the density $,, than $5. Several methods, varying 


from exact to simple estimations, are available to describe this effect. 


For simplicity the donor cell or upwind method is chosen here. This procedure looks at 
the direction of the convective flow and assumes that the value of a quantity at the interface 
is its value in the upstream control volume. Mathematically, the donor cell method evaluates 


$ at the west interface by 





(0.5 + 090, + (0.5 — 006, 


|| 


| 4.5 
Where: 


a = 0.5 if the flow 15 from W to O 


a= —(.5 if the flow is from Q tọ W 


& is the upwinding factor. It is evaluated using the time n value of the velocity whose sign 


indicates the direction of flow. If density is not defined on the interface, the procedure is 


followed again using velocity as the flow reference. Combining the results above, Eq. 4.1 is 


transformed to: 


ntl n 
Сеа = source HP Varia), (0.5 + a, oy" +(0.5-0,)00" ) 


Ну а) (0.5 + 0900 - (0.5 - o,)6z" ) 


Круя) (0.5 + a, Jê; * 4 (0.5 — 0,)05" ) 


n n 


Hp И. у ((0.5 + о,)%)” | T (0.5 ni OL, Ox" ) 4.8 
Rearranging and making the following substitutions 


M_=(p'v'a), 

A. =(0.5+0,)M, 
A. 2 — (0.5 — o, M, 
A. «(0.5 + о )М, 


ay = (0.5 > о,)М, 
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reduces Eq. 4.6 to 


n +1 А 
ve) = source -9p (A, +A +A, +A +M -M +M „— М) 


+A +А фу +A Y нА фт. 4.7 


This is the discretized form of the general equation. However, to put it in a more tractable 
form the operator H’ and the coefficient Ag are defined and substituted into the equation. 


The H’ operator is also used in several of the other discrete conservation equations. 


Ay=-(4,+4, +4, +4, +M_-M,+M,-M) 
H'Q) -A S; A Qw FAR FAGS 


° SH -4, =source +H’( ) ESAE 4.8 





Eq. 4.8 is discretized with the exception of the source term which is specific to the bal- 
ance being studied. Note, that the superscript * has been used in lieu of nel. Since the 
momentum and energy balances compute only an approximation of the n+1 value, the single 


asterisk designates it as the first guess. 


4.3.2 The Momentum Balance and Source Terms 


The discrete form of the momentum balances can be obtained by substituting the 


appropriate velocity for $ and transforming the source terms. 
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Fig. 4.2 shows a control volume associated with radial velocity and locates where the 
principal variables are defined. The pressure force acting as a source of radial momentum is 


discretized by estimating the derivative and multiplying by the area. 
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Figure 4.2 
The Radial Momentum Balance 


Control Volume and Variables 
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dP A 
pressure force = – 2а = (Р-Р; )а 4.9 


The current time step vatues of pressure are used since the PISO algorithm solves for veloc- 
ity using the momentum conservation equations before advancing the pressure variable. The 


mean gas phase area is used as discussed in section 3.3.2. 


The other source term in the momentum equations is the friction force. In the frits and 
the fuel particles this is characterized by the Ergun relation, Eq. 3.11. Rearranging this to 


calculate the pressure force yields 


150uv.(1-e% 1.75pvolvl(i-e 
Ara “ай 9 pA e a шс 4.10 
D € D, £ 


Since the friction force is a function of velocity squared, it must be linearized. The same 


technique used in the flux term is applied. Using Fig. 4.2 


n n n n n 2 n 2 | 
|у | 2 NC25(v4uw + Vane + Vase t+ Vasw)) + (У,о) 4.12 


This maximizes the implicitness of the relation. As a result, the left hand side of Eq. 4.8 


becomes 





р"бу 150u (1-8) | l?75plvol (1-2) || | p'óV - Т 
лонац КОРЕЕ СРИИ 
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The friction term for the plenums is treated similarly with the Re calculated based on the time 


n velocity. 


The discrete form of the radial momentum conservation equation used in the PISO 


solution algorithm is: 


De: ae. epit p vV,o9V 
pon es p vo tH Qu D er аи 4.14 


Eq. 4.14 has been arranged to facilitate solution of the resulting system of simultaneous equa- 
tions when each of the radial velocity control volumes is considered. The left hand side can 
be factored into a matrix of coefficients multiplied by a vector containing the radial velocities 
for each of the control volumes. The right hand side 1s a constant vector, since all values are 


known for time n. The axial momentum is handled analogousiy. 


4.3.3 The Energy Balance and Source Terms 


Heat transferred from the solid phase and the change in the pressure volume product 
are the source terms in the energy balance for the control volume shown in Fig 4.3. The 
energy transferred from the fuel particles couples the solutions of the solid and the gas phase 


models. 


The heat source term 1s simply: 
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Figure 4.3 
The Energy and Mass Balance 


Control Volume and Variables 


70 









> at ра 
Я 
o 
. К В © 
mv BLZ - 
77 і 09; fara se, 


ú 
ello 4 пу le IOA "n 


h @tugiq " 


onaligi m ан мены | и. 


asidaisav опа беи | ео. 









B 





ú E 


Q =; ҺАУ, (7; = (0 4. l 5 


The problem is whether to use the time n or n+1 values for the temperatures. A heat source 
term based on time n would be the easiest to handle. However, this introduces the stability 
problems associated with explicit differencing. These are magnified by the small gas volume 
and heat capacity. The authors of the THERMIT code (R-1) developed a fully implicit 


method for treating the heat transfer from the fuel which is used in this model. 


In order to treat both the solid and gas temperatures implicitly, the procedure consists 


of three steps which combine the results of the solid and gas energy conservation equations. 


First, the model assumes T7 equals T7 '' and solves the solid phase energy balance for 


T" and OT; *'/9T¿*' where T” is the first approximation to the new solid temperature. 


Second with T" substituted for T? * ' the fluid phase enthalpy balance is solved for the 


enthalpy at n-1. To do this, Eq. 4.15 is modified to account for the fact that T" is not the final 
value of Tz *'. This results in Eq. 4.16a which is transformed to the conservation variable 
used in the rest of the equation by substituting 5.16b. Eq. 5.16c shows the final form of the 


heat source term in the gas phase model. 
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* n+ oT. T , 
o navare anaw [Е lag -To 4.16a 


Qna 
Я и АЛ" 
ке се 4.16b 
P P 
Q МА Е T we A" - м f Р Дои h" x ie 
= ) ——+— T+ 7 | —— l| س‎ ! 
i s C» 6: = ý Ор: Cp Cp i 


Lastly, now that the true value of Tg *' is known, T' is updated to T; *' using Eq. 4.17. 


n+] 


$ "+1 n 
ELLAS -TO 4.17 
ОНА я 


Chapter 5 presents the fuel particle model and how to evaluate T" and the temperature 


derivative. The fuel temperature and heat transfer coefficient are further modified and 


replaced by T and U when the exact definitions of the solid phase model are applied. 


The pressure volume term is less complicated. The discretized form is 


OV(P 1. 
y- ( )о 


4.18 
dr Of 


Here (P./.). is the estimate of the pressure change between n and n+1 calculated by the first 


corrector in the PISO method. This value is available, since the energy balance is used as the 


second predictor. 
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Combining the enthalpy source terms with the general discretized balance equation 


results in the following discrete energy conservation equation: 


бур" UV_A О ВЕЕ д" апо 
ы Зад ръси » = UAVS ЕЕ ee 
of Cp dl ç J Cp ОВ" Кр 


+Н (А) 4.19 








n+] 


p is the first approximation of p" *' which is calculated in the first corrector step and Ay is 


the particle surface area per unit volume of the packed bed. 


4.4 THE MASS BALANCE 


The mass balance is implemented by using the control volumes centered on pressure 
and the other state variables. This in conjunction with the staggered grid, which conveniently 
defines the velocities on the faces of the control volume as shown in Fig. 4.3, makes the con- 


servation of mass flow calculation straight forward. 
Since no sources exist within the control volume, only the mass flows and the time 
derivative terms of the continuity equation require transformation. The techniques used on 


the similar terms in the momentum balances yield Eq. 4.20 which is the form of the mass 


balance used in the first corrector. 


127 bv +(pv,a), -(pv,a) + (ру, а = (ру, а) =0 4.20 
t 


To 





+ 


: п+ 1. : 
Note that the second estimate of v^" is used. No requirement for a mass balance was 


imposed on the first predictor. As a result, the v values do not satisfy this relation. 


4.5 THE PRESSURE EQUATION 


The pressure equation is the core of the PISO method. It is the relation that forces both 
mass and momentum to be conserved simultaneously. The derivation of the equation out- 


lined below is adopted from references (I-} ,[-2). 


The continuity equation, 4.20, for the control volume and variables shown in Fig. 4.3 
forms the foundation of the pressure equation. To this several substitutions developed from 
the momentum equations are made. Focusing on the east face axial velocity for a moment, 
the first momentum predictor and corrector relations take the form of Eq. 4.21 and 4.22 


respectively. Note that the corrector form is explicit. 


OV Ao * * OV 

ч 4.21‏ + 5 — ? == " کے ی ل 
E " } ve=HWw.)-(APa)+tp Vue s‏ 

oV Ao . ж ж e бу 

а в р и 4.22 
| St p" } УЕ Н (Var) ( a), TP VE бї 


Subtracting Eq. 4.21 from 4.22 and solving for p'v;; one obtains: 
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p vi -K (AP “a » oe K^ (AP"a), +) v 4.23 


ta 
Of р” 


АР” =Р.-Р, 





Multiplying by the area then gives the mass flow rate through the east face of the control vol- 
ume. Substituting this and the similar expressions derived from the velocities at the other 


faces yields: 


oV ж n P rt 
-=r (P 0 AO TID tD, ы?) 
+О ДРЕ - РЕ)+0 (Ру - Ру) 
+ D,(Py — Py) * DP; = Р;) =(р"у,) = (р"у,), +(р"у,), — (p^v,), 


D. z(K a^), 4.24 


The last step relates the density difference to the pressure using the equation of state, 


P =pRT, to obtain the following: 





4.25 


This combined with Eq 4.24 and some algebra gives the final form of the first corrector pres- 


sure equation. 
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The pressure equation used in the second corrector is derived in the same fashion. Using Eq. 
4.22 and the second corrector for momentum, Eq 4.27, gives Eq. 4.28 which is the second 


corrector pressure equation. 


oV Ao «жо wee š жж “с non OV 
eS УЕ =Н (v) - (AP a). +p VaE 5; 4.27 
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Where: 
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The pressure equation roots in the momentum and mass conservation equations force 
the pressure and velocity fields to satisfy the two conditions simultaneously, thereby giving 


one the option of using a noniterative solution procedure. 


4.6 THE SOLUTION ALGORITHM 


Sections 4.3 through 4.5 developed the equations required for the PISO solution 
scheme outlined in section 4.2.2. The following steps provide the solution sequence used for 


the gas phase in the fuel element model. 


First Predictor 
1. Solve the radial momentum balance, Eq. 4.14, for the radial velocity field, v, . 


2. Solve axial momentum balance, Eq. 4.14, for the axial velocity field, v,. 


First Corrector 
3. Solve the pressure equation, 4.26, for the pressure increment, P' — P". 
4. Calculate the P field by adding the pressure increment to P". 
5. Use the equation of state as modeled in the H2EOS program to compute p` 
given P` and T". 
6. Compute the new radial velocities, v; , using Eq. 4.23. 


7. Calculate the new axial velocity field, v; , using an equation similar to 4.23. 
Second Predictor 


8. Obtain A from the enthalpy balance. 


9. Use the H2EOS program to get T;; given P^ and A". 
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Second Corrector 
10. Solve the second pressure equation, 4.28, for P" — P”. 
11. Calculate P". This is now P"*' 
12. Compute p" using the equation of state given P" and Tz. This isp" ^ *. 
13. Calculate the radial velocity, v,  , using a relation equivalent to Eq. 4.27. 
This is the time n+1 value. 
14. Calculate the axial velocity, v, , using a relation similar to Eq. 4.27. This is 


the n+1 time step value. 


The equations solved in steps 1, 2, 3,8 and 10 are actually sets of M times N simulta- 
neous linear equations, where M is the number of radial nodes and N is the number of axial 


nodes. These steps are most easily accomplished by treating the equations in matrix form. 


[Matrix of Coeff. ][ Variable Vector] = [Constant Vector] 


A wide variety of methods exist to solve this problem. For large numbers of nodes iterative 
techniques are probably most efficient. Issa recommends the use of ADI or Stone’s Strongly 
Implicit Procedure (I-2). For smaller matrices direct solution methods such as LU decompo- 


sition can be used. 


4.7 BOUNDARY CONDITIONS 


Once the equations and physical data are available, problem definition is completed by 
adding the appropriate boundary conditions. Figures 4.4 a, b and c show the control volume 


arrays for typical momentum and pressure equation solutions. 
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Figure 4.4a 
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Control Volume Array 
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Figure 4.4b 
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Figure 4.4C 


Pressure Control Volume Array 





The user determines the inlet pressure and temperature (PBI and TBI) and the outlet 
pressure (PBO) when the transient is defined. The other boundary conditions are imposed by 
the fuel element geometry. One constraint common to all equations is that the velocity per- 
pendicular to a solid surface is 0 on the surface. Thus, the radial velocities that would be 
defined on the outer surface of the inlet plenum (B in Fig. 4.4a) are known to be O and the 
points are not included in the solution matrix. Symmetry also requires this condition on the 


fuel element center-line. 


A constraint similar to the radial case exists for the axial momentum at the west and 
east ends, Fig 4.4b. However, the arrangement of the fuel element and control volumes 
impose some complications. The first occurs in the region of the hot frit. The design pre- 
vents axial flow in the frit, since it is constructed from a solid piece of metal drilled to allow 
radial flow. The second occurs along the west end. The inlet and outlet port control volumes 
are separated by a series of control volumes which contain half solid boundary and half par- 
ticle bed. The axial velocities are defined on the solid surface and are therefore 0. The prob- 
lem is how to maintain an array of control volumes with a solution matrix that retains the 
appropriate coupling and mass balances, yet still gives correct solutions for the boundary 
values. 

Two methods exist to handle this problem. One is to eliminate the control volumes for 
known boundary values. The resulting nonrectangular array complicates the programing if 
the user still has the option to define the number of axial and radial nodes. It also risks 
decoupling the outlet plenum from the rest of the fuel element. The other method, the one 
used here, includes the boundary volumes necessary to achieve a rectangular control volume 
array and forces the solution to the correct value by controlling the source terms in the 


boundary control volumes. Using an artificially high friction term effectively reduces the 
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included boundary velocities to 0. This has the advantage that anisotropic friction character- 
istics can be used to adjust the axial velocity with minimal effect on the radial velocity. As a 


result radial flow is still allowed in the west end volumes and the hot frit. 


One other boundary value issue that affects the two pressure equations is the evaluation 
of the pressure gradient across the boundary control volumes. The model assumes that there 
is no difference between the pressure at the boundary surface and the pressure in the center of 
the control volume. If the pressure at the center of the control volume changes, so does the 
pressure at the solid surface. This gives a zero pressure gradient at the boundary. The pres- 
sure gradient across the opposite side of the control volume can still vary in accordance with 
the variable pressure on either side of the interface. This is consistent with the lumped 


parameter approach. Figure 4.5 illustrates the principle. 





NN 
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Figure 4.5 


The Boundary Pressure Gradient 


48 SUMMARY 


Discretization and linearization of the control volume balance equations results in sets 
of simultaneous equations which can be solved for the pressure, velocity, temperature and 
density in each contro! volume. The PISO method is used to accomplish the numerical solu- 
tion. It uses a 2 stage predictor corrector technique which advances the variables from one 
time step to the next by calculating a sertes if successively refined approximations. Its 


semi-implicit, noniterative structure provides ease of solution and numerical stability. 
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CHAPTER 5 
THE SOLID PHASE MODEL 


5.1 OBJECTIVE 


Chapters 3 and 4 present the gas phase part of the model. Equally critical to the suc- 
cess Of the fuel element code is the modeling of the fuel particles. They are the source of the 
heat and contain the majority of the stored energy. This chapter presents the development of 


the solid phase model equations and the solution method. 


3.2 THE MODEL 


The solid phase is simpler than the gas phase, since no motion is present. As a result, 
predicting the solid phase temperatures requires only a heat balance. The energy conserva- 
tion model must account for three processes, heat deposition, storage and dissipation. The 
first two can be described on the scale of a single fuel particle. The third represents the fuel 


particle's interaction with surrounding particles and the hydrogen coolant. 


The lumped parameter, control volume approach used in the gas phase model is also 
used here. Because this method assumes all particles in a control volume are the same, the 
single particle processes can be analyzed on an individual basis, then rescaled to account for 
the total solid mass in the control volume. The lumped parameter approach is also conve- 
nient for the dissipation mechanisms of conduction and radiation. These processes are 
dependent on a temperature gradient. Since no gradient exists within a control volume, they 
can be considered as occurring between control volumes in the overall energy balance. Eq. 


5.1 is the general form of the energy conservation equation for a control volume. 
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k — effectivity conductivity 


q — heat source 


5.2.1 The Fuel Particle 


The temperature distribution within a fuel particle is a function of the heat deposition 
rate, heat removal rate and fuel particle dimensions and materials. The most accurate model 
would divide a fuel particle into at least four nodes (based on composition), then create and 
solve a set of finite difference equations. This process would have to be repeated for one 
particle in each control volume, making the model extremely cumbersome. To simplify the 
particle model, the single node approach is adopted. This reduces the fuel particle to a single 


composite control volume which retains all of the properties of the original. 


Reference M-2 derives the single node model for a cylindrical fuel pin. The same tech- 
nique is used here except that spherical geometry has been substituted. The single node anal- 
ysis homogenizes the fuel particle by assuming a steady state temperature distribution and 
constant material properties in each region of the fuel particle. The resulting particle heat 


balance takes the form of Eq. 5.2 
карен К, ру 
Шр T-2q'-U(T-T,) 9:2 


Where: 


т, = particle mass per unit outside surface area (kg/m‘) 
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Cp = average particle specific heat (J/kg К) 


T = effective fuel particle temperature (K) 
47 - heat generation per unit outside surface area (W/m) 


U = effective overall heat transfer coefficient (W/m? K) 


The second term on the right side represents the heat transfer to the gas phase. The model 
code uses the following equations to evaluate the terms in the single node relation. These are 
easily derived using the procedure of reference M-2. The number subscripts refer to regions 


identified in Fig 5.1 
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Figure 5,1 
The Fuel Particle 
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T is, in effect, defined to be the temperature that satisfies the requirements of the single node 


heat balance. It is not a true average or associated with any particular point in the fuel par- 
ticle. The determination of the heat transfer coefficient, h, is discussed in chapter 3. These 
relations assume that the average temperature of a region in the fuel particle can be 


approximated by the arithmetic average of the regions inner and outer interface temperatures. 


The material properties of the fuel particle components vary with temperature. This 
model assumes fixed values for the product of the density and specific heat. Thermal con- 
ductivities are allowed to vary with temperature and are evaluated using the correlations of 
reference D-2. Even with these correlations, the values of the physical parameters are only 
estimates. For example, estimates of the carbide fuel conductivity vary by as much as 30% 


between references. 
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The single node analysis also pennits easy estimation of the heat flux time constant, T. 


If power is increased linearly, Tau is the time delay between achieving a given power level 
and having its equivalent steady state flux on the outside surface of the particle. This 
assumes that the coolant temperature remains constant (M-2). Fig 5.2 diagrams the relation- 
ship. As such, it is a measure of the time for the energy to move from the fuel kernel to the 


surface of the particle and into the gas phase. 





1 = 5.4 


In this case, t is a function of coolant flow and temperature and varies from 30 to 105 ms. 


9 
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5.2.2 Solid Phase Conduction and Radiation 





Heat transfer between adjacent solid particles occurs via several mechanisms. Conduc- 
tion can take place through particle to particle contact points or through the gas trapped 
between the particles near the points of contact. If temperatures are sufficiently high, 
radiative heat transfer also contributes (Y-1). The effective conductivity of the solid phase is 


the sum of these components. 


The dominant mode of conduction is through the fluid near the point of contact. 
Because the fluid is in the boundary layer, its conduction properties are relatively insensitive 
to the gas flow rate except at very high Reynolds numbers (Y-1). The packed bed element 
Reynolds numbers are low because of the high fuel particle surface area to bed volume ratio. 
Direct conduction via the points of contact is only significant at low pressures. As a result, 
the conductivity of the solid phase is highly dependent on the fluid and is therefore a function 
of the gas and solid conductivities. Kunii and Smith proposed the following relation to deter- 


mine the effective conductivity. It is the one incorporated in the fuel element code. 


ВА -=). 
EE un "o | 5.5 


ks 


p+; 
Where: 


k° = effective bed thermal conductivity 
ko = gas thermal conductivity 


€ = bed void fraction 
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В = packing parameter, assumed to be one for a randomly packed 


bed 


ф = a measure of the effective thickness of the liquid film between 


particles near the contact point 


$ is calculated using the following equation: 





(= sin” 0 p 
== o so. = ^^” ШЕТ. ER из ар 26 
NESR NCO OE СОЗ) 
Ks 
K = " 
l 
Ө=— 


n is the number of points of contact with other particles and represents the number of possi- 
ble paths for conduction. Because the number of points of contact is not known exactly ina 
randomly packed bed, ó is calculated for the most open packing, 6, , and a maximum close 
packing, €, . The values are averaged, weighting them by the true bed void fraction as shown 


in Eq. 5.7. 


€ — .0260 


9 —(6,— 95) 0.216 t, 5.7 
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The temperature required for significant radiative effects is a function of particle size. 
For 1 mm particles thermal radiation is important at temperatures greater than 700 K and 
1,800 K for .] mm diameter particles (S-4). With anticipated fuel temperatures of 2,000 K 


radiative heat transfer is significant and is included in the model. 


Schotte includes two paths in his model of the radiation effects. One is the radiation 
across the void space to an adjacent particle, k^, and the other is the radiative transfer 
between particles in series with conduction through the particles. The applicable equations 


for SI units are: 


3 


| Т 

k^ 20.229(EM)£D, 0° 5.8a 
1-£ 

= 5.8b 
ou 


F 


k, 1s the effective radiation conductivity of the bed and EM is the emissivity of the fuel par- 


ticles. The fuel element code uses the properties of the outer layer, zirconium carbide, for the 
conduction parameters because it is the material in contact with the other particles and its 
surface characteristics determine the emissivity. The code uses an emissivity of .78. The 
radiation component of conduction is a function of the surface temperature as well. The 
model program approximates this with the effective solid temperature, 7. The actual surface 
temperature is less than T but the size and spherical shape of the particles minimize the AT 


within the particle, mitigating the error of the model. 
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Combining the conduction and radiation conductivities yields the overall effective con- 


ductivity of the solid phase, Kora. 
ku =k +k, 5.9 


5.2.3 The Solid Phase Solution Method 


The solid phase energy balance is solved by applying a discretized form of the energy 
balance to an array of control volumes, just as is done with the gas phase model. The control 


volumes are coincident with their gas phase counterparts. 


Assimilating the discussion of the previous sections into the energy balance requires 
rescaling the single node analysis of the heat deposition, storage and transfer to the gas phase 
from an isolated particle to a control volume basis. The terms of Eq. 5.2 were derived on a 
per unit particle surface area. The rescaling is accomplished by multiplying both sides of the 
product of the control volume size, V,,, and particle surface area per unit volume, Ay. The 


units for each term in Eq. 5.10 then become Watts. 


_ d — zu E 
M,C V Av T =q" V Ay- UV A(T Т) 5.10 


The area per volume for spheres is a function of the particle diameter and the void fraction. 


mo E) 
an 


Av 5.11 
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The last step in assembling the continuous form of the control volume energy conserva- 
tion equation is to add the conduction terms which account for the heat transferred through 


the control volume boundaries. This yields: 


E "uo E E kt 
nm, С» „Аус T = q V.A, = OV AMT ==: To) ey > oat УТ, > 1 2 


Note that the area used in the conduction term is the area of the control volume face, not the 


surface area of the particles. 


Equation 5.12 must be discretized to make it compatible with numerical solution tech- 
niques. Since the control volumes are the same as those of the gas phase, the solid tempera- 
tures are defined at the same points as the gas state variables. Solid temperatures are not 
defined in the plenum areas. No solid exists in these contro} volumes and conduction to the 


moderator block is not considered. 


The discretization procedure is the same used in chapter 4. The last term on the right 
hand side is the only new term which has not been discretized earlier. The temperature gra- 
dient across an interface is approximated by the difference in the temperatures on either side 
of the interface divided by the distance between the temperature points. The problem now is 


how to evaluate the conductivity at the interface. 


Conductivity is a function of temperature, void fraction and material. Therefore, sig- 
nificant differences may exist on opposing sides of the interface. The first impulse is to use 
the arithmetic average of the two values. However, this has been shown to give erroneous 


results (P-2). The harmonic mean is a better choice and has several physical arguments in its 


a 





favor. Patankar equates the heat flow on either side of the interface to derive the relation for 
the harmonic mean (P-2). Eq. 5.13 shows the relation for the east face of a control volume. f 
is the ratio of the distance between the east interface and T, to the distance between the cen- 


ters of the two control volumes. Similar relations can be derived for die other faces. 


ЖШ 
=| —+—— 
Кы), Kr ko 5 . | 3 


- 


The source term representing the heat transfer from the solid to the coolant Is treated 
implicitly as discussed in chapter 4. Therefore the n+! values of the temperatures are used. 
In accordance with the THERMIT implicit method, 77. is used as an estimate of the T7; '. A 


correction is made later in the solution process. 


Substituting the discrete forms of the terms into Eq. 5.12 results in the following form 


of the heat balance: 


5 TO y Ж ое; (Т жа 
т СУ А EN ШЕП та лайна рог. osa pur e SER 
| | ЕСІ) Ж р 
к шлык A ` 
total d Sr s total or ; 
+q” V „Ay - ОУ АТ" Т) 5.14 


This can be rearranged into a form similar to the gas phase balances with these substi- 


tutions: 


98 














kaqa = А. Koad м А, К,а E A, Kora E A. 
2 da ór J ór J. 


ИАА В 





ef 


n+} —n +) 


\=A7T, ЗАТ, tA Ry +A 


п +1 


H “( P 5.15 


Equation 5.16 is the final form of the solid phase heat balance used in the fuel element model 


code. 
-(B -Ao - UV „A To +H (T ) 2 4 V, A,— UV, AVTz - BT, 5.16 


Setting up this equation and evaluating the coefficients for each control volume creates 
a system of simultaneous linear equations. The model code can now solve these for T with 
the same matrix solution techniques used for the gas phase. The values obtained for the solid 


temperature are only a first estimate. 


The THERMIT method for solving the coupling of the gas and solid phases also 


° : : 4 mn +i . . š 
requires the evaluation of the temperature derivative, OT" /dTZ*' . Continuing with the 


л + 


assumptions that Tc and T equal T2 *! and T | respectively, the code differentiates Eq. 5.16 


to get: 


Ae +1 m 


` 5 
ОТ (B-A,-UV,A,) 
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Before completing the solution for the solid phase temperature, the gas phase energy 
conservation equation must be solved for the true value of 72°' . When this is done, all of 


the terms in Eq. 5.18 are known and T" can be corrected to get T í 


=n +! =e or’ n+ " 
Г = "оре Чо = G) 53.18 


This process ıs equivalent to saying 
Interphase Heat Transfer = СУ АТ" Ti) 5.19 


in both the gas and solid phase equations. 


5.3 SUMMARY 


An energy balance is applied to the fuel particles to model heat generation, storage, 
and transfer. Transfer includes convection to the coolant conduction to surrounding particles 
and radiation to adjacent material. To simplify the model a steady state temperature distribu- 
tion within the fuel particle is assumed. This allows the state of the solid phase to repre- 


sented by a single temperature in each control volume. 
The coolant and fuel particle models are coupled by the interphase heat transfer term. 


In order to preserve numerical stability, the model treats both gas and solid temperatures 


implicitly using the technique developed for the THERMIT code. 
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The model of the packed bed fuel element is now complete. 
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CHAPTER 6 
RESULTS AND CONCLUSIONS 


6.1 OBJECTIVE 


Chapters 2 through 5 present the details of a mathematical model of a packed fuei ele- 
ment. It is based on the momentum, energy and mass conservation equations. These rela- 
tions are encoded in a micro computer Fortran program described in appendix B. As a partial 
test of the model's validity, the code is applied to a series of steady state and transient 


problems. This chapter presents the results of these test cases. 


6.2 STEADY STATE 


6.2.1 Problem Statement and Methodology 


Six steady state problems test the model's ability to reach a reasonable solution under 
various conditions. Changing the peak power level provides the simplest means of simulat- 
ing the wide range of hydraulic and thermal conditions that the fuel element would experi- 
ence during pulsed high power operation. Cases with peak power densities of 2.1 and 1.0 
GW/m' are examined as well as a zero power problem. The other boundary conditions are 


held constant to permit more meaningful comparisons. 


In addition to power density the operator must also specify the following: 
1. The physical dimensions of the fuel element 
2. The control volume dimensions and arrangement 


3. The initial values of all variables 
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4. The coolant inlet temperature and pressure 
5. The initial and final outlet pressure 


6. The initial and final power density 


All of the test problems are conducted using a rectangular array of control volumes 
arranged with 5 equal sized nodes axially and 10 nodes radially. The radial nodes vary in 
size to match the fuel element components. One node was assigned to each plenum and frit 
and 6 nodes were assigned to the fuel particle bed. The control volumes and variables are 
identified by two indexes, I and J for the radial and axial positions respectively. Figure 6.1 
shows the indexing scheme. The origin of the grid is the pressure control volume at the inlet 
port (1,1). Note that the velocity control volumes cross component boundaries, because of 
the nature of the staggered grid. For example, radial velocity (2,1) is defined on the interface 


between the inlet plenum and the cold frit. 


The sign convention is positive for axial velocities moving from low index to high. 
Radial velocities are positive for flow from the outside of the fuel element toward the center. 
Thus, increasing I and J represent the positive directions. Under normal flow circumstances 
then, axial velocity is positive in the inlet plenum and negative in the outlet plenum. Radial 


velocity is positive in the fuel bed under these circumstances. 


For steady state cases the initial conditions are arbitrary. However, they should repre- 
sent an equilibrium condition to ensure initial numerical stability. The test problems are 
started with steady flow, zero power and an element pressure drop of 50 kPa. All 


temperatures are 300 K and inlet pressure is 2 MPa. The final outlet pressure is 1.915 MPa. 


for all cases. 
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VR(i+1,J) 


VZ(1,J) VZ(1,J+1) 





VR(1,J) 


3 State Variables (pressure, temperature) 
VR Radial Velocity 


VZ Axial Velocity 


Figure 6.1 


Staggered Grid 


Indexing Scheme 
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The code takes the peak power density and applies a chopped cosine weighting to sim- 
ulate the axial power profile. In this case, the factors were .89, .97, 1.0, .97, and .89 for axial 


nodes | through 5 respectively. 


The cold frit thickness varies with position to match the axial flow profile and the 
power distribution. The thicknesses are based on the PIPE element design drawings. Values 


of 2.26, 1.93, 1.80, 1.78 and 1.91 mm are used for nodes | through 5 respectively. 


The problem is completely defined without having to specify the mass flow rate. It 


varies as the conditions in the fuel element change. 


The steady state solution is actually achieved by conducting a transient for a long 
enough period of time that the values of the variables become essentially non-varying. The 
transient 1s initiated by making a step change in the outlet pressure and power density to the 
final values. This essentially simulates instantaneously raising power and opening the fuel 


element outlet valve. 


Steady state is identified by monitoring the overall system mass and heat balances. All 
problems are run for 3.75 sec. which insures that the energy transferred to the coolant is 
within 5% of the energy deposited by the simulated fission. The mass balances agree within 


1% or less. 


The test cases also provide a test of the numerical stability of the solution technique. 


As expected, the semi-implicit differencing used in the discretization limits the size of the 
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time step that can be used. Instabilities can be initiated by raising the power or increasing the 
fuel element pressure drop. A time step of .25 ms works well for the power densities and 


pressures used in the test problems. 


6.2.2 Hydraulic Results 


The principle hydraulic issues address the subjects of model complexity, | or 2 dimen- 
sions, and flow distribution prediction. Figure 6.2 shows the calculated steady state values of 
the variables for the 2.1 GW/m' case (Refer to Fig. 6.1 to meld the arrays). The zeros on the 
edges of the velocity arrays represent velocities defined on the fuel element boundaries. The 
hot frit shows zero axial velocities due to the modification of the axial friction term to simu- 
late the frit's drilled single piece construction. The small diameter of the particles in the cold 


frit also result in high friction and negligible axial flow. 
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Steady State Fuel Element 
Power Density = 2.1 GW/m’, Inlet Femp.— 300 К, 
Inlet Press. = 2.0 MPa, Outlet Press. = 1.915 MPa 


Superficial Radial Velocities (m/sec) 


11| .0000 .0000 000 BODO .0000 

101 6.6141 6.0407 6.0424 544790 3.3305 
25397 3.715. 5.7147 5.4678 5.0468 
8| 5.1095 4.7442 4.7188 4.5146 4.1770 
7| 4.0823 3 [837 3.8043 3.6624 3.3877 
6| 3.1475 2.9651 2m 2.0951 2.6506 
5| 2.3003 2.1792 2.1830 2.0985 1.9494 
4| 1.5059 1.4474 1.4315 1.3704 1225 
3j .5821 ‚3724 3428 Sil 4685 
2| .5109 5004 4703 4392 995 
ll .0000 0000 .0000 .0000 .0000 





Superíicial Axial Velocity (m/sec) 


10 


= NU E CA G. і оо О 





-223.3224 


0000 
.0000 
.0000 
ОООО 
.0000 
.0000 
.0000 
.0000 

6.4169 


-173.2407 


.0000 
-1.4026 
-1.5122 
- 1.6399 
OT 
-1.9537 
-2.0629 

-.0017 
5.0037 


-127.9471 


.0000 
-.9275 
- 9552 

-1.0684 
SEVIK 
-1.3024 
-1.4397 
-.0010 
3.6201 


-82.7846 


.0000 
-.5075 
-. 3323 
-. 5699 
-.0236 
-.7041 
-.8257 
-.0005 
2.3196 


Pressure (MPa) 


10 


= М ә Б л С ~) CO 





ЏЈ 


1.92157 
1.92359 
1.92610 
92711 
1.92788 
1.92845 
1.92884 
1.92913 
1.96743 
2.00001 


1.93258 
1.93461 
1.93716 
1.93822 
1.93904 
1.93966 
1.94010 
1.94045 
107 77 
2.00002 


2 


1.93969 
1.94179 
1.9444] 
1.94544 
1.94623 
1.94681 
1.94721 
1.94749 
1.97619 
2.00003 


3 


107 


1.94380 
1.94584 
1.94835 
1.94931 
1.95004 
1.95056 
1.95090 
1.95112 
1.97798 
2.00003 


4 


Figure 6.2 Steady State Results - 2.1 GW/m’ 
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In the fuel bed axial velocities are significant and vary from 17 to 360% of the radial 
velocity, with the maximum just inside the cold frit. The associated flow is toward the ported 
end of the fuel element (Jzi). This implies that the high pressure drop of the cold frit forces 
the coolant to distribute itself over the entire surface of the cold frit. Once through the frit, 
the coolant redistributes as it follows the many parallel paths through the fuel bed and outlet 
plenum to the exit port. The axial flow increases from the closed end to the ported end of the 


element. 


Figures 6.3a,b and c show the mass flow profiles for the three cases, comparing the 
axial distribution at the cold frit and the hot frit. At all power levels, the flow shifts to the 
ported end of the element after passing through the cold frit. This is consistent with the axial 


velocity field. The redistribution increases with power. 
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Flow Distribution 
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Axial Position (node) 


Figure 6.3a 
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Flow Distribution 


10 GW/m"3 


Mass Flow (%) 
30 


25 


NO 
о 


= 





Í —F —— з‏ — ——— کیہ م 


........ 
..... 24... 
......... 


... 
D00 
... 





Axial Position (node) 
Figure 6.3b 
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Axial Position (node) 


Figure 6.3C 
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The flow distribution through the cold frit is also a function of power. Figure 6.4a 
illustrates the shift from a centrally peaked flow profile at zero power to one peaked in the 
first axial node at 2.1 GW/m’. The hot frit exhibits similar behavior only it is more pro- 


nounced, Fig. 6.4b. 


The sources of the pressure losses in the element also change with power density. Fig- 
ures 6.5 and 6.6 plot the changes in the pressure drop along the length of the outlet plenum 
and across the cold frit as a function of power. At zero power, 74 to 82 kPa of the 85 kPa 
total element pressure drop occurs across the cold frit. However, at 2.1 GW/m' the cold frit 
represents only 48 to 71 kPa of the pressure drop. The biggest change is in the closed end of 
the element (J=5). The cold frit changes result from increases in the pressure losses in the 
fuel bed, hot frit and outlet plenum as the hydrogen heats up. This is expected, since the 
lower density of the warmed coolant forces higher velocities to maintain the mass flow rate. 
Therefore, frictional losses, modeled as a function of velocity and velocity squared, increase 


rapidly. 
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Figure 6.4b 
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Figure 6.6 
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Figure 6.7 
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Axial Position (node) 





The explanation of the hydraulic behavior of the model lends itself to an electrical anal- 
ogy. The fuel element can be viewed as an inlet and outlet port connected by an infinite 
number of parallel flow paths with different resistances, similar to parallel resistors in an 
electrical circuit. Just as the current distribution adjusts to give the same voltage drop across 
all of the resistors so does the mass flow distribution adjust to give the same pressure drop 
across all possible flow paths. The complication in the fuel element is that resistances are a 
function of power and therefore variable. As the hydrogen heats up, the relative resistance of 
the flow paths changes forcing a change in the flow profile. These changes favor increasing 
the percentage of flow in the more direct paths despite the greater thickness of the cold frit. 
The higher axial velocities in the fuel particle bed at high power represent increased flow tn 
the paths that bypass the outlet plenum. However, it should be noted that the vast majority of 
the flow is still collected in the outlet plenum and that the total shift is less than 10% in the 


cases studied. 


6.2.3 Thermal Results 


Thermal issues examined in the test problems include the temperature profile and what 
heat transfer processes should be included in the model. 

The temperature distribution along the hot frit shows significant variation, contrary to 
the design objective of approximately equal temperature gasses entering the outlet plenum, 
Fig.6.9. This is the result of the flow redistribution causing a mismatch in the mass flow and 


axial power profiles. 
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Figure 6.10 shows the radial coolant and solid phase temperature distributions in the 
J=3 node for a peak power density of 2.1 GW/m?. Looking at the J=1, 3 and 5 fuel control 


volumes next to the hot frit, the model predicts the following heat balance results: 


Control Volume (8,1) (8,3) (8,5) 
Heat Generation (W) 19849 22275 19849 
Power transferred to the 15685 14870 12203 


gas phase (W) 


The heat transferred to the gas phase is calculated using the overall heat transfer coeffi- 
cient and the difference in the fuel and coolant temperatures. If the heat balance is main- 
tained and steady state exists, then the difference between the heat generation and the power 
generation must represent the energy conducted or radiated to the surrounding fuel particles 
and the hot frit. The elevated temperatures in the cold frit are also the result of conduction 


from the fuel particles. 


The gas film temperature difference is a maximum in the first fuel particle control vol- 
ume inside the cold frit (I=3). The zero temperature difference in the cold frit is indicative of 
its excellent heat transfer characteristics, especially the high surface area to volume ratio 
(700,000 compared to 7,200 m’/m’ in the fuel). The smallest film temperature drop occurs 
next to the hot frit and is the result of conduction which reduces the power transfer required 


to the coolant. 
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The heat transfer tume constant for the single node model also varies with temperature. 
At 300 K it is approximately 100 ms in the fuel particles, while at 1935 K it is 30 ms. Look- 
ing at the formula for the time constant, this reflects an increase in the overall heat transfer 


coefficient, U. 


To further check the effect of conduction and radiation, a case without conduction 
between control volumes was run. Figure 6.11 shows the difference in the solid phase tem- 
peratures with and without conduction included in the model for the J=3 node. The fuel tem- 
perature next to the hot frit changed relatively little. However, fuel temperatures in the rest 
of the bed decreased significantly. The coolant temperature follows the fuel. Near the hot 
frit, the lack of conduction forces higher gas film temperature drops to accommodate the 
larger direct heat transfer from the fuel to the hydrogen. This was confirmed by the steady 


state heat balance where the heat transferred to the coolant equaled the power generated. 


The no conduction case also displayed a small increase in the mass flow rate from 36.7 


to 37.8 g/sec. The change in flow resulted in a 25 K decrease in the outlet temperature 
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6.3 TRANSIENT RESULTS 


A rapid up power transient is conducted to test the codes ability to track rapidly chang- 
ing conditions. In this case, steady state flow is established then the peak power is ramped 


from 0.0 to 2.0 GW/m' over a period of one second. 


Figures 6.12 through 6.16 show the time variation of the variables. The results are con- 
sistent with the trends observed in the comparison of the steady state results at various power 


levels. 


The power plot, Fig. 6.12, also provides some insight into the characteristics of the 
model fuel element. During the first second of the transient while power is increasing at a 
constant rate, the time constant of the system can be estimated from the time difference of the 
curves, since they are almost parallel to each other (See section 5.2). In this case the time 
constant is probably between .6 and 1.0 seconds which agrees with the times required to 
reach steady state. This is only an approximation, because the heat transfer parameters 


included in the time constant are themselves functions of temperature. 
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6.4 CONCLUSIONS 


While the data used in the test problems are representative of the Sandia PIPE experi- 
ment fuel element, the numerical results may not be directly applicable due to approxima- 
tions made where PIPE data are not readily available. However, from a more general 


perspective several basic ideas are demonstrated in the model results. 


Two dimensional modeling is required to reflect accurately the axial flow in the fuel 
particle bed which is present to some degree in all cases. Mass flow is a function of tempera- 
ture for a given pressure drop across the fuel element. It decreases as power increases due to 


the increased flow resistance of the hotter, lower density hydrogen coolant. 


As power increases, the pressure losses across the hot frit and in the outlet plenum also 
increase. Thus, if the total AP remains unchanged, then the pressure loss in the cold frit must 
decrease to compensate. This reduces the total mass flow rate and drives a partial redistribu- 
tion of the coolant flow to the shorter flow paths which see less of the change in the outlet 
plenum pressures. The redistribution is manifested in increased axial flow through the fuel 
bed and an increased percentage of the flow through the thicker sections of the cold frit near 


the inlet port. 


The cold frit flow distribution design is specific to the conditions assumed in the design 
process. If power pressure or flow rate are changed the cold frit will not provide the desired 
flow profile. Higher power favors the more direct flow paths from the inlet to the outlet port. 


This potentially could place the maximum flow in the area of minimum power. 
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The cold frit does not act in isolation. Therefore, its design should be approached from 
a system perspective that looks at the total fuel element pressure drop and the contributions 
from each component in each flow path. Once the cold frit design is fixed, it is the variation 
of the pressure losses in the other components with changing conditions that determines the 


final flow distribution. 


The thermal model demonstrates the necessity to consider conduction and radiation 


processes when predicting fuel element behavior. 
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СНАРТЕК 7 
SUMMARY AND CONCLUSIONS 


7.1 OBJECTIVE AND METHODOLOGY 


Considerable work is in progress to design and develop new space based reactors. The 
designs by necessity must be lightweight and small with high power densities. One reactor 
being studied at the Brookhaven and Sandia National Laboratories uses a hydrogen cooled 
packed bed fuel element. The MIT interest in this system is the development of a closed loop 


digital controller for reactor power. 


In order to control the reactor properly, the controller must evaluate the reactivity 
effects of the hydrogen coolant and fuel element temperature distribution. Therefore, in 
addition to a good neutronic model, a real time thermal hydraulic analysis capability is 
required. The purpose of this investigation is to perform the initial steps in the building of a 
thermal hydraulic module for the controller. Specifically, a detailed model of the fuel ele- 
ment thermal hydraulic behavior is presented. It forms the reference case against which sim- 


pler real time models can be compared. 


The first step in the modeling process is to look at the details of the actual fuel element. 
This model is based on the element constructed for the Pulsed Irradiation of a Particle-bed 
Element (PIPE) experiment being conducted at the Sandıa National Laboratory. The cross- 
section is shown in Fig. 7.1. The fuel is installed in the annular space between two cylindri- 


cal concentrically mounted porous retention elements (frits). Coolant is directed around the 
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outside of the outer frit then turns and flows radially through the cold frit, fuel and inner frit 
to the collection plenum at the center of the element. The fuel consists of 500 um diameter 


particles with a UC, core and pyrographite and zirconium carbide coatings. 


The form of the model is determined by the identifying the features and processes that 
affect the fuel element behavior. The flow portion must accurately represent the velocity and 
pressure of the coolant as a function of position. Thus, manifold distribution effects and fric- 
tion processes are important. Coolant and fuel particle temperature are also required. Asa 


result, intra and interphase heat transfer processes must be considered. 


The transport mechanisms can be related by fundamental principles of momentum, 
energy and mass conservation. Applying these with a lumped parameter approach to an 


array of control volumes forms the basis of the fuel element model. 


The unifying structure is provided by using a staggered grid, Fig 7.2, to link the control 
volumes. The staggered grid is based on placing points where velocity is defined between 
the points where pressure is defined. In this arrangement, the control volumes can be 
matched to the components of the fuel element. The fuel element model assumes cylindrical 


symmetry, allowing the problem to be reduced to two dimensions (radial and axial). 
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The final step in constructing the model is the discretization of the control volume con- 


servation equations to make them compatible with digital computer solution techniques. 


The gas and solid phases are modeled independently except for the interphase heat 
transfer term which couples the two models. The coupling is addressed in the solution pro- 


CESS. 


7.2 THE GAS PHASE 
7.2.1 The Model 


The gaseous hydrogen coolant is highly compressible. As a result, sumultaneous 
momentum, mass and energy balances are required to define the state of the coolant. The 
momentum and energy: balances are identical in form and can be described using the general 


variable 6 . 
d 
q f= È sources + Zit 7.1 


ф represents the superficial radial and axial velocities in the momentum balances and 


the specific enthalpy in the energy balance. 
The momentum source terms include the pressure drop across the control volume and 


the frictional losses. The pressure drop is easily evaluated, since the staggered grid defines 


pressure on the faces of the momentum control volume. 
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Two models are used for the friction term. The frits and the packed bed are character- 
ized by the Ergun relation which treats flow in the bed as flow through tubes with irregular 


cross-sections (E-1). 
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From Eq. 7.2 it can be seen that frictional losses are a function of velocity and velocity 
squared. v, is the component of velocity in the direction of the momentum balance and | v| is 


the magnitude of the superficial velocity. 


Frictional losses in the plenums are described using a pipe flow relation that assumes 


they are proportional to a friction factor, f. 
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L = length D = equivalent pipe diameter 
Eq. 7.3b is a modified form of the McAdams relation to account for surface roughness. 


The cold and hot frits and their associated plenums serve as dividing and combining 
manifolds respectively. Bajura and other investigators found that direct application of the 
momentum conservation equations to manifolds did not accurately represent experimental 
behavior because they failed to account for the axial momentum contained in the portion of 


the stream which turns in the control volume and leaves it flowing radially (B-4,B-5 and 
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B-6). They found that adding two constants to the momentum balance equation which char- 
acterize the velocity profile of the axial stream and the axial momentum of the turning stream 
significantly improved the accuracy of his model. A similar analysis was applied to the 
outlet plenum. The momentum equations are summarized in Table 7.1 where the constants 


have been combined to form the turning coefficient, Ө. 


The energy balance contains two source terms. The first is a pressure change term. 
The second enthalpy source is the heat transferred from the fuel particles to the coolant. It is 
a function of the temperature difference across the gas film, the overall heat transfer coeffi- 
cient and the surface area per unit volume ratio. The heat transfer coefficient is evaluated 
using the Colbum j factor empirical correlation and the single node analysis of the fuel 


particles (B-3). The final fonn of the enthalpy balance is shown in Table 7.1. 


Table 7.1 


The Gas Phase Conservation Equations' 
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The mass balance, also shown in Table 7.1, uses only the flux term, since there is no 


source of mass in the control volume. 


7.2.2 The Solution 


Before the variables can be advanced in time using the balance equations, the continu- 
ous form of the equations must be discretized. The fuel element model equations are devel- 
oped using a semi-implicit five point difference method. A donor cell, upwind, system is 
used to evaluate the scalar variables at the control volume interfaces. In the discrete form the 
balances are implicit using the time n+! value of the balance variable. However, the coeffi- 
cients (such as the mass flow in the flux term) are evaluated using the time n values. This 
introduces a semi-implicit nature to the overall equation. The discretization also serves to 
linearize the relations. In cases where squared terms appear, they are approximated by using 


the product of the time n+1 and time n values keeping the equation linear in the n+1 variable. 


The Pressure Implicit with Splitting of Operators (PISO) method forms the basis of the 
solution algorithm (I-1, [-2). As used in this model, it is a two stage predictor corrector tech- 
nique that advances the variables in time by calculating up to three successively more refined 
estimates of the n+1 value of the variables. Because it is noniterative and has reasonably 
good stability, it is relatively fast and requires less computing effort than competing methods. 
Each step of the solution sequence applies the appropriate balance equation to each control 
volume to generate a system of simultaneous linear equations which are solved using stan- 


dard matrix techniques. 


The first predictor solves the momentum balances for the velocity field. Then taking 


advantage of the staggered grid, which defines the velocities on the faces of a pressure cen- 
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tered control volume, the momentum and mass balances are combined to form a pressure 
equation which yields new velocities and pressures that simultaneously satisfy both balances. 
The second predictor ıs the energy balance which contains the coupling to the solid phase 
solution. The second corrector applies the pressure equation again with the updated tempera- 
tures. At the end of the process three estimates of the velocity field, two of the pressure and 


one of the temperature have been made. 


7.3 THE SOLID PHASE 


7.3.1 The Model 


The solid phase is abstracted to mathematical form by using the energy balance. The 
balance includes three basic processes, heat generation and storage, energy transfer to the 


coolant and heat transfer to the surrounding solid material. 


The first two processes are initially modeled on a single particle basis. To simplify the 
model, the single node technique of reference M-1 is used. This assumes a steady state tem- 
perature distribution in the multilayer fuel particle and then homogenizes it to calculate effec- 
tive overall properties which are representative of the original particle. The temperature 
distribution in the particle can then be described by a single temperature which satisfies the 
heat balance equation. Since the lumped parameter approach is being used, only one particle 
needs to be modeled in each control volume. The results can then be rescaled to include all 


of the particles in the control volume. 


Heat transfer to the surrounding particles is modeled on a control volume basis. The 


driving force is the temperature difference between adjacent control volumes. Two mecha- 
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nisms contribute to the transfer. The conduction process їп a packed bed actually occurs 
through the gas layer near the points of contact. The Kunii and Smith model is used to 
evaluate the control volume conductivities (K-3). In addition to thermal conduction, radia- 
tive mechanisms are significant when particle temperatures exceed 1000 K. This is evaluated 
using the relation proposed by Schotte (S-4). The effective conductivity of the packed bed is 
the sum of the thermal and radiative components. Because of differences between the con- 
trol volumes the harmonic mean of the conductivities is used to estimate the effective con- 


ductivity at the interface. 


The solid phase energy balance then becomes: 
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This is discretized in the same fashion as the gas phase relations. 


7.3.2 The Coupling of the Gas and Solid Phase Solutions 


The last part of the model is the coupling of the gas and solid phases. The heat trans- 
ferred between phases is the only term common to both models. The relatively small volume 
and heat capacity of the hydrogen require that this term be treated implicitly for numerical 


stability. This is accomplished using the technique developed by Reed and Stewart for the 


144 





THERMIT code (R-1). It first estimates the new fuel temperature assuming the old value of 
the gas temperature. The method then solves the coolant energy balance using the estimated 
solid temperature plus a correction for the original assumption, based on d7;/dT, which is 

calculated in the solid phase energy balance. Finally, the estimate of the solid phase temper- 


ature is corrected. 


7.4 MODEL VERIFICATION AND RESULTS 


The model was encoded in a micro computer based Fortran program and applied to 
several steady state and transient problems to test its validity. Required input information 
includes the physical dimensions of the element, initial values of the variables, the nodaliza- 
tion scheme and boundary conditions. The boundary conditions can be completely specified 
by the inlet pressure, inlet temperature, outlet pressure and peak power density, all of which 
are easily measured quantities. Transients can be modeled by varying the outlet pressure 


and/or the power density. 


The steady state cases were identical except for the power density which was varied 
from zero to 2.1 GW/m’. The hydraulic results revealed significant axial flow in the fuel par- 
ticle bed which resulted in a redistribution of the coolant mass flow between the frits as 


shown in Fig. 7.3. This effect is more pronounced as power increases. 


Increasing power also changes the distribution of flow along the inlet plenum and cold 
frit. This is illustrated in Fig 7.4. This second redistribution effect was also noted when the 
model was altered to prevent axial flow in the fuel bed. Power also affected the mass flow 
rate. Flow decreased when power increased for a given pressure drop across the fuel ele- 


ment. 
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These effects are the direct result of heating the coolant, reducing its density, thereby 
increasing its velocity and causing a realignment of the pressure field. In the zero power 
case, 85 to 95% of the pressure drop occurs across the cold frit. When power is raised to 2.1 
GW/m', only 62% of the pressure drop occurs across the cold frit (at the axial center of the 
element). The frictional losses increase due to the higher velocities, especially in the hot frit 
and outlet plenum. Looking at the fuel element as a network of parallel paths between the 
inlet and outlet ports (the paths must have a common pressure drop), the higher velocities in 
the bed and outlet plenum increase the pressure drop in the paths through these elements. In 
order to satisfy the total pressure condition, which remains constant, flow must decrease in 
the high loss paths. This is reflected by an increase in the flow through the more direct paths 
between the ports and an increase in the axial flow in the bed that bypasses the outlet plenum. 


The total redistribution was less than 1046 in the cases examined. 


Although the thickness of the cold frit is varied to force a particular mass flow distribu- 
tion, the flow can be significantly altered by the outlet plenum and hot frit. The cold frit does 
not act in isolation. The thickness must therefore be chosen to obtain a specified flow pattern 


at power in concert with the design of the other fuel element components. 


The thermal results predict significant temperature variation along the hot frit, contrary 
to the design objective of approximately equal temperatures. This is the direct result of the 
mass flow redistribution producing a mismatch in the power and flow profiles. The tempera- 


ture distribution is shown in Fig 7.5. 
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The test cases also indicate that conduction and radiation are important mechanisms for 
heat transfer especially at the higher fuel temperatures near the hot frit. In these locations up 
to 30% of the energy generated in the fuel particles is removed from the control volume by 


transfer to adjacent solid material in the fuel element. 


A transient problem, which raised peak power from zero to 2.0 GW/m' in one second, 
confirmed the trends observed in the comparisons of the steady state results. The system 


time constant was estimated at .6 to 1.0 seconds. 


7.5 RECOMMENDATIONS FOR FURTHER INVESTIGATION 


The model developed in this investigation represents a first step in the design and con- 
struction of a thermal hydraulic module of a controller for reactor power. It is the reference 
case against which other simpler and faster versions can be compared. Several refinements 


would greatly increase its utility. 


First, the model predictions should be compared to experimental results. Several 
model parameters such as the turning coefficients in the momentum equations may require 
adjustment, since the current values, based on the literature, were obtained under much dif- 


ferent conditions. 


Second, the reference case should be corrected and expanded. The incorrect use of the 
superficial velocity instead of the interstitial velocity in the momentum balance time deriv- 
ative and flux terms should be corrected. Preliminary calculations based on steady state 


results show that the error in the momentum flux term is approximately 3% of the pressure 
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source term in the momentum equations shown in Table 7.1. The effect of the error on tran- 
sient calculations has not been evaluated. Within the fuel element the effects of conduction 
to the moderator block and component expansion due to heating have not been addressed. 
The latter could have a significant effect on the hot frit and fuel bed pressure losses. Once 
the fuel element analysis is complete, it should be broadened to include the entire core and 


reactor system. 


The third area for investigation starts the process of designing a real time controller 
module. A series of model simplifications should be developed and tested against the refer- 
ence case to determine which provide the required prediction accuracy. Areas that should be 


looked at are: 


1. The numerical scheme. Is the PISO algorithm really the most efficient numerical 
technique or would the technique of one of the commercial codes be better? 

2. The nodalization. Is a fine node analysis really required or can the fuel element be 
adequately represented with just a few control volumes? Can one fuel element be modeled in 
detail and the results generalized to the other elements in the core? 

3. Model basis. Can the model be simplified by using another approach to the solution 
of the basic conservation equations? For example, could the momentum integral method 


used in the MINET code be applied (M-3,V-2,V-3)? 


Packed bed reactors have tremendous potential because of their high power density and 


inherently safe design. It is hoped that the thermal hydraulic model presented here will con- 


tribute to their development. 


151 





APPENDIX A 
FUEL ELEMENT MODEL 
PROGRAM DESIGN AND OPERATION 


AT OBJECTIVE 


HTWOCOOL is a PC compatible program that performs the calculations associated 
with the model developed in chapters 3, 4 and 5. This appendix presents the operational 


details required to use the code for further analysis. 


A.2 PROGRAM DESCRIPTION 


The HTWOCOOL code advances the fuel element pressure, velocity and temperature 
fields in time using initial values of the variables and a set of boundary conditions. The code 


is capable of tracking transients or calculating steady state solutions. 


The code is written in Fortran 77 using the Microsoft 4.10 compiler. As such, it is 
intended for use On a personal computer. However, the test problems conducted to date indi- 
cate that the time step must be restricted to values between .25 and .5 ms. This small time 
step requires exceedingly long computing tumes for PCs (48 hrs or more on an AT 
compatible to get to steady state). Therefore, conversion to a mini computer may be more 


practical than using a PC 


The program’s design is modular with a shell program to assemble data and sequence 


the subroutines. The following modules comprise HTWOCOOL: 
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Shell Program 
MAIN: Assembles data, sequences the subroutines, keeps track of the time step, 
calculates transient parameters. 
Solid Phase Subroutines 
ONCNST: Calculates constants for the single node fuel particle model. 
SOLID: Performs the solid phase heat balance. 
PISO Subroutines 
AXLMOM: Performs the axial momentum balance. 
RADMOM: Performs the radial momentum balance. 
PRESUR: Solves the first stage pressure equation. 
ENTHLP: Calculates the coolant energy balance and corrects the initial estimate 
of solid phase temperature. 
PRESII: Solves the second stage pressure equation. 
Support Subroutines 
H2EOS: Evaluates the parahydrogen physical constants. 


LUDCMP and LUBKSB: Matrix equation solution routines. 


To advance the variables one time step MAIN establishes matrices for each variable 
with its time n value in each control volume. It then passes the appropriate data to each of 
the subroutines which calculate the successive approximations to the time n+1 values. Each 
of these approximations is maintained in its own array. When the calculation for the time 
step is complete and the results recorded, MAIN substitutes the n+1 values into the time n 


matrices to set up the calculation of the next time step and advances the time counters. 
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The PISO and Solid phase subroutines function in a similar manner. They take the 
appropriate conservation equation for each of the control volumes and evaluate the coeffi- 
cients of the balance variable using the time n values received from the shell program. This 


creates a set of simultaneous linear equations in the following form. 


[Coefficient Matrix][Variable Vector]=[Constant Vector] 


The subroutine then calls LUDCMP and LUBKSB including the matrix and constant 
vector as arguments. These routines perform an L-U decomposition and back substitution 
respectively to solve for the new value of the variable. This result is passed back to the call- 
ing subroutine in the form of a modified constant vector. The PISO or Solid subroutine then 
converts the vector to a two dimensional matrix of the updated values and updates other 
variables as appropriate. For example, PRESUR calculates the new pressure field then 
advances the velocity and density fields. Finally, the results and control are passed back to 


the shell program. 


The updated variables are identified by the number of Ts in the name. For example VZ 
and VZO are the original axial velocity field. VZST, VZSTT and VZSTTT are the first, sec- 
ond and third updates of the velocity. The third estimate is the one carried forward to the 


next time step. Similar naming schemes have been set up for the other variables. 
The nodalization scheme is defined by the user in the input data. The code does require 


that one radial node be assigned to the each plenum and frit. Any deviation from this will not 


work. However, there is no restriction on the dimensions of these preassigned nodes. 
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The sign convention used in the development of the code is positive radially moving 
from the outside toward the center-line of the fuel element. Axially velocity is positive mov- 
ing away from the inlet end of the element. The north side of a control volume is the side 


facing the center-line. The other faces line up according to the compass. 


The H2EOS module was developed by the National Bureau of Standards and modified 
by the Sandia National Laboratory (R-2). The LUDCMP and LUBKSB routines were taken 


from reference P-5. 


A.3 PROGRAM OPERATION 


The program is simple to run. It is not interactive. The only operator action required is 


the assembly of the data set. All operational choices are entered via the data set. 


Power and outlet pressure transients can be run as well as steady state calculations. 
Transients are either in the form of a step change or ramp and are initiated immediately or 
after a time delay. In the case of ramped change the final value is specified and the time over 


which the transient is to occur (TOPT or TOQT). The applicable equations are 





Time 
Power Density = A +B 
PTinte 
Outlet Pressure = де ра De 


where: 
QTime =Time — Delay time 


B=Final Power Density 
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Note that the overall power density is the actual variable used in the transient simu- 
lation. Once MAIN determines the power density for the time step, it applies a chopped 
cosine function to represent the axial power distribution. Power density is assumed to be 
constant radially. Changes to the power distribution function would require program modifi- 


cations and recompiling. 


Once the data file is prepared, the program is executed by entering the word MAIN. A 
MAKE file is also included should program modifications be required. This file is a Micro- 
soft utility that contains the compile and link commands. After modifying the source code of 
the any of the routines, simply enter "make main" (quotes omitted) and the utility will 


recompile the updated modules and relink the program giving a new executable file. 


A.4 INPUT-OUTPUT 


The input is list directed and the majority of the output is formatted. The list directed 
input leaves the operator some freedom to place decimal points and vary the size of the 
entries. The values are entered sequentially, in the order the program reads them. Each entry 
Is separated by a comma. There is no restriction on the number of values per line except that 
the program automatically skips to the next line of data when a new READ statement is 


executed. 


The SI system of units is used exclusively. All lengths are entered in meters, masses in 


kilograms, pressures in Pascals and temperature in Kelvin. 


The Program looks for the input data in a file called PBR.DAT. This file contains the 


following information: 


156 





a. Physical dimension of the cell, fuel and control volumes. 

b. User defined boundary conditions. 

c. Time step size, number of time steps to be executed. 

d. Initial conditions, TGO, TSO, PO, VRO, VZO matrices. | 


e. Transient control parameters. 


To define the control volumes the program requires the axial length (DZ) and the a vec- 
tor of radii. The axial dimension is assumed to be the same for all nodes. The grid is vari- 
able in the radial dimension. While considerable latitude in the fuel element and node 
configuration is allowed, some of the geometry has been built into the program. The nodes 
for the coolant state variables align with the components of the fuel element. The velocity 
nodes, however, may lie across component boundaries (see Fig. 2.6 The Staggered Grid). Up 


to ten nodes radially and 10 nodes axially can be defined for a total of 100 control volumes. 


The fuel element is pictured to be lying horizontally so that in an I, J array of control 
volumes the values of I and J represent the radial and axial positions respectively. I follows 
flow, increasing as the position moves inward. J starts at the end of the fuel element with the 


inlet and outlet ports and follows flow in the inlet plenum. 


The values of the radius vector define the centers and surfaces of the control volumes. 
R(1) is the outside surface of the inlet plenum and R(2*M-1) is the center-line of the fuel 
element, where M is the number of radial nodes. The radial position of the center of the Ith 
node is designated R(2*I) and the north and south faces are R(2*I+1) and R(2*I-1) respec- 


tively. 


157 





The program assumes a minimum of 4 nodes radially. Any additional nodes are added 


to the fuel region of the element. The fixed control volume assignment assignments are: 


_I _ Location 
1 Inlet Plenum 
2 Cold Frit 
M-1 Hot Frit 


M Outlet Plenum 


Table A.1 gives the order of data entry and a line breakdown. 


Output of the velocity, pressure, density, enthalpy and temperature variables is for- 
matted and prints as an array. The staggered grid shows the relationship of the variables. 
Note that the program as written inverts the order of the points showing the ceolant entering 
at the top of the fuel element rather than at the bottom as shown in chapters 4, 5 and 6. The 


indexes printed next to the array are correct. 


Also included in the output are list directed records of the time step, power level, over- 
all heat balance and overall mass balance which facilitate following a transient. Refer to the 
MAIN program listing to determine exactly which ones are used. Some of the output is 


directed to the screen. 
The formatted output is accomplished by the RESULT subroutine. This may be called 


at any time during the program, after every time step or at the end of the transient. Program 


modifications and recompiling are required to change it. 
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The current program is configured to direct the results to an output file and the screen. 


However, Fortran does allow the option of going directly to the printer. The required OPEN 


statement is included in the program but not activated. 


Line 


10 


ІЙ 


Table A.1 

Input Data Format 
Data Elements (in required sequence) 
Number of radial nodes (M), number of axial nodes (N), inlet coolant 
temperature, inlet pressure, fuel particle diameter, cold frit particle diam- 
eter, axial node length, time step size, number of time steps to advance 
variables 
Void fraction (M values starting with [=1) 


Cold frit thickness (N values starting with J=1) 


Radii to control volume centers and interfaces (Starting with outside of 
inlet plenum and working inward to 2*I+1; this may take several lines) 


Initial value of the radial velocity (The order should be all the J=i values 
from l=1 to M+1, then all the values for J=2 etc to J=N. This may be 
more than one line.) 


Initial value of the axial velocities (The same order as the radial veloci- 
ties except I=1 to M and J=1 to N+1) 


Initial value of pressure (The same order as the radial velocities except 
that I=1 to M and Jz1 to M) 


Initial solid temperature (Use the same entry scheme as pressure) 
Initial coolant temperature field (Use the same sequence as the pressure) 


Pressure transient parameters: initial outlet pressure, final outlet pres- 
sure, transient duration, time delay for initiation of the transient 


Power transient parameters: initial power density, final power density, 
transient duration, time delay for initiation of the transient 
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APPENDIX B 
PROGRAM LISTING 


The following program listing contains the source code for the original segments of the 
program. The code for the H2EOS, LUDCMP, LUBKSB modules is contained in references 


R-2 and P-5. 
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(ОТО ОСЕ 


E RE (AAA 


PROGRAM TO EVALUATE THE HYDROGEN COOLANT TEMPERATURE PRESSURE AND 
DENSITY IN A PACKED BED REACTOR 


INITIALIZATION AND DATA INPUT 


COMMON/GEOTRY/R (45) , EP (20) , DZ, THKCF (30) ,DV(20),DP,DPCF,M,N,DT 
COMMON /BOUCON/PBI, PBO, TBI, RHOBI, HBI 
FEWT DT D2, VRI11, 10). v2 (10, 12), P (1067 T0 TG 10, 10) TS (10/10) 
+ RHO(10,10),RHST(10,10),VZST(10,11),VRST(11,10),VRSTT (11,10), 
House (10) 11) ,BRAP (11,10), B2AF (10, 11) 
REAL VRSTTT (11,10), VZSTTT (10,11), TIME,A,B,C, TAU, PWRDEN 
REAL PST(10,10) ,O(10, 10), PIO(10, 10) 
REAL BVR(11,10),ARF(11,10),AWR(11,10),AER(11,10),ASR(11, 10) 
REAL ANR(11,10) 
BER 28,2 (10, 11,,AzZPi(10, 117, AW2(10,12) 20210107 a82110, 11) 
REAL ANZ(10,11) 
REAL H(10,10),R,EP,PBI,PBO,RHOBI,TBI,PI 
ERATEAOZ(1I0,11) Б2(10 11) АОВ(11,10) БЕ(11,10), ТИКСЕ 
BEAR HST(10,10).TGST(10, 10), PSTT(10, 10) , RBSTT (19, 10) 
ПЕВНО СТ ТО) У20(10,119,/Р0(10, 10), Т50(10 10)  ТСО(10/10) 
REAL RHOO(10,10),HO(10,10),HTI, HTO, NN, KK 
REAL PBOI, PBOF, TOPT, PDLAY, PTIME, MAVE 
REAL QTIME, TOOT, QDLAY 
REAL MA(10),CPBAR (10) ‚UT(10,10),FBAR(10,10),TSST(10,10) 
REAL UBAR(10,10),AV(10) 
REAL OIN(10, 10), TSSTT (10, 10) , ANN(10, 10), KZC (10, 10) 
WMRECER M, N UNITS, Т, J, IT, TSTOP, СТ 
PARAMETER (PI=3.14159, UNITS=1) 
NOTE: TO CHANGE THE GRID SIZE THE TYPE STATEMENTS MUST BE MODIFIED 
IN ADDITION TO THE DATA FILE. 


ASSIGN PRINTER, INPUT DATA FILE AND OUTPUT DATA FILE 

OPEN (3, FILE-'LPTI1') 

OPEN (4, FILE=’C: \LANG\FORBIN\PROG\PBR. DAT’ ) 
OPEN (4, FILE=’ PBR.DAT’ ) 
OPEN (3, FILE=’ PBRRES.DAT’ ) 

IMPORT GEOMETRY DATA AND CONSTANTS 
READ (4,*}M, N, TSI,PB1,DP,DPCF,Dz, DT, TSTOP 
READ (4,*) (EP(I), I=1,M) 

READ (4,*) (THKCF (I), I=1,N) 
READ (4,*) (R(I), 1=1,2*М+1) 

TMPORT INITIAL GUESS VALUES FOR VARIABLES 
READ (4,*) ((VRO(I,J), I=1,M+1), J=1,N) 
READ (4,*) ((VZO(I,J), I=1,M), J=1,N+1) 
RESO (4, *) (BOL IS), 11M), ет N) 

READ (4,*) ((TSO(I,J), I=1,M), J=1,N) 
READ (4,*) ((TGO(1,J), I=1,M), J=1,N) 

IMPORT TRANSIENT PARAMETERS 
READ (4,*) PBOI, PBOF, TOPT, PDLAY 
READ (4,*) A,B, TOQT,QDLAY 

TNITIZLIZE TIME 
TIME=0. 
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RHOBI=PTDENS (PBI, TBI,UNITS) 
HBI=PTENTH(PBI, TBI,UNITS) 
NN=REAL (N) 

CT=0 

DO 565 TT=1,TSTOP 


CONSTRUCT THE INITIAL MATRICIES FOR RHOO, HO AND CONTROL 
INITIALIZE MATRICES USED FOR CALCULATING COEFFICIENTS 


DO 200 I-1,M 
DO 100 J-1,N 
TE (1 EOS 1) TCO(1, 9) TBI 
TE (тт Ст Ту сото 10 
RHOO (I, J) -PTDENS (PO(I,J), TGO(I,J) , UNITS) 
HO (I,J) =PTENTH (PO(I, J), TGO(I, J), UNITS) 
Б(ғ 2)-БО(1 7) 
TS (I,J)2TSO(I,J) 
TG CISJ)STGOlI 9) 
RHO (I, J) =RHOO (I, J) 
H(I, J) =HO(I,J) 
CONTINUE 
DV (I) =PI* (R(2*I-1) **2-R(2*I+1) **2) *DZ*EP (I) 
CONTINUE 
DO 202 I-1,M*1 
DO 201 J-1,N 
VR (1, J) =VRO (I, J) 
CONTINUE 
CONTINUE 
DO 204 I=1,M 
DO 203 J=1,N+1 
V2 (1,J)=V20(1,9) 
CONTINUE 
CONTINUE 


CALCULATE TRANSIENT PARAMETERS 


TES (TIME GE: BDLAY) THEN 
PTIME=TIME-PDLAY 
ELSE 
РТІМЕ-0. 
ENDIF 
TE A TIME GE. QODLAY) THEN 
QTIME=TIME-QDLAY 
ELSE 
QTIME=0. 
ENDIF 
WRITE (6, *)” TIME’ , TIME 
WRITE(6,')^PTIME',PTIME,'OTIME',OTIME 
ТЕСІРТІМЕ „СТ. ТОРТ) СОТО 206 
TE (TOP T EQ О.) TREN 


PBO=PBOF 
GOTO 206 
ENDIF 


PBO=- (PTIME/TOPT) * (PBOI-PBOF) +PBOI 
WEITE(6, *) 7 PRO’ PBO 

IF (QTIME .GT. TOQT) GOTO 207 
PWRDEN=A+B*QTIME 
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207 WRITE (6, *)’PWR DENSITY’ , PWRDEN 
DO 220 L=3,M-2 
DO 210 K=1,N 


KK=REAL (K) 
Q(L,K) =PWRDEN* (DV (L) /EP (L) ) * (COS (2*PI-(1~((2*KR~1}/NN) ) 
+ *.5886)) 
210 CONTINUE 
220 CONTINUE 


С PERFORM FIRST UPDATE OF SOLID PHASE TEMPERATURE 
CALL ONCNST (TSO, MA, CPBAR, UT, FBAR, KZC) 
CALL SOLID (MA, CPBAR, UT, FBAR, TSO, TGO, PO, RHO,Q, VRO, VZO, KZC, 
+ TSST, UBAR, AV, ANN) 
E CALL RESULT(VR,VZ,PO,RHO,TG,H, TSST) 
E PERFORM GAS PHASE UPDATE 
CALL RADMOM (VR, VRO, VZ, RHO, RHOO, P, TG, BRAF, VRST, BVR, ARF, AWR, 


+ AER, ASR, ANR) 

300 GALL, AXEMOM(VRE, VZ, VZO, RHO, RHOO. P, TG, BZAF, VZST, BVZ, AZF, AWZ, 
+ AEZ, ASZ, ANZ) 

E CATT RESULT (VEST, VIST ©. RHO. Te. HTS) 

500 CALL PRESUR (VRST,VZST, RHO, RHOO, RHST, P,PO, PST, P10, TG, BZAF, 
x BRAF, VRSTT, VZSTT) 

G GALÍ RESULT (VRSIT, VZST.. PSt RHST, 0G, H, ro) 


Ghul ENTHLP(VZSTI, VRSTT, RHO, RHOO, RHST, PST, P10,H, H0,UBAR, TSST, 
Wane GO, ANN, HST, TGST, AV,OIN, TSSTT) 
6 CALE RESULT (VRSTT, VZSIT, PST, RHST, T1GST, HST, TS STT) 
CALL PRESII (VRST, V2ST, VRSTT, VZSTT, 2, PST, PO, TEST, RECO, 
+ RHST, VRSTTT, VZSTTT, PSTT, RHSTT, BVZ, BVR, ARF, AZF, AWR, 
: AER, ASR, ANR, AWZ, AEZ, ASZ, ANZ, RHO) 
C | HEAT BALANCE CHECK 
QCONV-O. 
HTI-0. 
DO 542 K-3,M-2 
DO 541 L-1,N 
HTI-O(K, L)+HTI 
QCONV=QCONV+QIN (K, L) 


541 CON TINUE 
542 CONTINUE 
HAVES(ABS(VZSTTT(M,1))*RHSTT(M,I)*PI*(R(2*M-1)**2) 
+ TABS(VZSTIT(1,1))*RHOBI*PI*(R(1)**2-R(3)**2))/72 
HTO=MAVE* (HO (M, 1) -HBI) 
E WRITE (3, *)’ HEAT BALANCE’, HTI, HTO, QCONV 


WRITE (6, *)’ HEAT BALANCE’, HTI, HTO, QCONV 
WRITE (6, *) QIN (M-2,1) , QIN (M-2, 3) , QIN (M-2, 5) 
WETIEBCO,*)O(M-2,1),0(M-2,3) ,Q(M^2, 9) 
БЕРІСТІ “ЕО: 20), "THEN 
IP TTE O TIME C HTO OCONV  VASITI(M I) RHSII(M ] IE LA (R(2%M=1)%%2 
+ J TGST(M 1) HTL, TSSTT(M=-2 SF  IOST(M=2;, 3) 
СТ=1 
ELSE 
CT=CT+1 
ENDIF 
ПЕЕ 6G ~)HST (1,3), PST (5,3), Pst. (575) 
560 CONTINUE 
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561 
563 
565 
600 
9001 


9002 


2003 
9004 


9005 


9006 
9007 


9008 


9020 


2009 
9010 


9011 


2012 
9013 


9014 


9015 
9016 


TIME=TIME+DT 
DO 563 E=1,M 
DO 561 K=1,N 


VRO (L, K) =VRSTTT (L, K) 
VZO (L,K)=VZSTTT(L,K) 
RHOO (L,K)=RHSTT (L,K) 
PO (L,K)=PSTT (L,K) 

HO (L,K)=HST (L,K) 
DOOCb K)eSTGST(L.K) 
TSO(L, K) =TSSTT (L, K) 


CONTINUE 
CONTINUE 
CONTINUE 
WRITE (3, 9001) 


FORMAT ('0',6X, 


WRITE (3, 9002) 


EORMAT('O*", IOX,"' 


, 6’) 


"RADIAL VELOCITY”) 


тко а тве, 8 3 тве е а" тве 157 TRE, 


DO 9004 I=1,M+1 
ПЕТ ТЕРС(З, ЗООЗУЛ/(УЕКӘТІТСІ 21-1) 


PORMAT (" 0", 6X, 


CONTINUE 
WRITE (3, 9005) 


FORMAT (’ 0’, 6X, 


WRITE (3, 9002) 
DO 9007 I=1,M 
WRITE (3, 9006) 
FORMAT ('0', 
CONTINUE 

WRITE (3,9008) 


FORMAT ('0',6X, 


WRITE (3, 9020) 


BORMAT(’ 0” , 10X, ’ 


, 3.) 
DO 9010 I=1,M 
WRITE (3, 9009) 


FORMAT (’0’, 6X, 


CONTINUE 
WRITE (3,9011) 


FORMAT (’ 0’ , 6X, 


WRITE (3, 9002) 
DO 9013 I=1,M 
WRITE (3, 9012) 


EORMAT (707 ,6X%, 


CONTINUE 
WRITE(3,9014) 


EORMAT(’0776X, 


WRITE (3, 9002) 
DO 9016 I=1,M 
NEITE(3, 2015) 


FORMAT (’ 0’ , 6X, 


CONTINUE 
WRITE (3, 9017) 


6X, 


12,3Х,5(Е8.4,3Х)) 


AXIAL. VELOCITY” ) 


IOUVZSTITTI SJ) JST NT) 
12,3X,6(F10.4,1X)) 


PRESSURE?) 


ое. 21 швед“ 3 TIES i 41. TAG; 


І, (РӘТТКТ,У),4-1,М) 
T2,3%,1P,61E11,5,1X%)) 


DENSITY?) 


I, (RHSTT(I,J),J=1,N) 
12, ЗХ, 5 (Еб.4, ЗХ) ) 


’COOLANT TEMPERATURE’) 


T.(TGEST(I,J},)=1,N) 
Ше ЗХ (Е. 4772 
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9017 


9018 
2019 


9021 


9022 
570 


590 


600 
9001 


9002 


9003 
9004 


9005 


9006 
9007 


9008 


9020 


9009 
9010 


9011 


9012 
9013 


FORMAT (’ 0’ , 6X, ' ENTHALPY’ ) 

WRITE (3, 9020) 

DO 9019 I=1,M 

WRITE (3,9018) I,(HST(1,J),J=1,N) 

FORMAT ('0',6X,12,3X,1P,5E10.4) 

CONTINUE 

WRITE (3,9021) 

FORMAT (’ 0’ , 6X,' SOLID TEMPERATURE’ ) 

WRITE (3, 9002) 

DO 9022 I=1,M 

WEITE (3,9015) 1, (TSSTT(1,9),90=1,N) 
CONTINUE 

STOP 

END 

SUBROUTINE RESULT (VR, VZ,P,RHO, TG, H, TSST) 
COMMON/GEOTRY/R (45), EP (20) , DZ, THKCF (30) , DV(20) , DP, DPCF, M, N, DT 
COMMON/BOUCON/PBI,PBO, TBI, RHOBI,HBI 
РЕТ Та) и (ТО ТТ) P (010.10) 95BHO (10410), TC (10. 10) (10 10) 
REAL R,EP,DZ, THKCF,DV,DP,DPCF,DT, TSST (10, 10) 
ВЕАГ, РВТ, РВО, ТВТ, ВНОВТ, НВТ 

INTEGER I,J,K,L,M,N 

WRITE (3,9001) 

FORMAT ('0',6X,' RADIAL VELOCITY’ ) 

WRITE (3,9002) 

HORMAT(^0' ,10X,* “1, ,TR6,' 2°, TRG,” 3’. тво,’ U TRG. 5°, TRG, 
, 6”) 

DO 9004 I=1,M+1 
WRITE(3, 9003) I, (VR(1,J),J=1,N) 
FORMAT ('0',6X,12,3X,5(F8.4,3X)) 
CONTINUE 

WRITE (3,9005) 

FORMAT ('0',6X,' AXIAL VELOCITY’) 
WRITE (3,9002) 

DO 9007 I=1,M 

WRITE (3,9006) I,(VZ(1,J),J=1,N+1) 
FORMAT ('0',6X,12,3X,6(F10.4,1X)) 
CONTINUE 

WRITE (3,9008) 

FORMAT ('0',6X,' PRESSURE') 
WRITE (3, 9020) 

FORMAT ('0',10X,' TRE 2',TR6,' 3' , TR6,' 4' , TR6, 
, 5”) 

DO 9010 I=1,M 

WRITE (3,9009) 1,(P(1,J),J=1,N) 
FORMAT ('0',6X,12,3X,1P,6(E11.5,1X)) 
CONTINUE 

WRITE (3, 9011) 
FORMAT (’0’, 6X, ‘DENSITY’ ) 
WRITE (3, 9002) 

DO 9013 I=1,M 

WRITE(3,9012) I, (RHO(I,J),J=1,N) 
FORMAT (’0’ , 6X, 12, 3X, 5(F6.4, 3X) ) 
CONTINUE 
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WRITE (3,9014) 

FORMAT ('0’,6X,'COOLANT TEMPERATURE’ ) 
WRITE (3, 9002) 

DO 9016 I=1,M 

WRITE (3,9015) I,(TG(1,J),J=1,N) 
FORMAT ('0',6X,12,3X,5(F6.1,4X)) 
CONTINUE 

WRITE (3, 9017) 

FORMAT ('0', 6X, ' ENTHALPY’ ) 
WRITE (3, 9020) 

DO 9019 I=1,M 

WRITE (3,9018) I, (H(I,J),J=1,N) 
FORMAT ('0',6X,12,3X,1P,5E10.4) 
CONTINUE 

WRITE (3, 9021) 

FORMAT ('0',6X,' SOLID TEMPERATUPE' ) 
WRITE (3, 9002) 

DO 9022 I=1,M 

WRite (3, 9015) 1, (TSST(I, J) , Jel, i) 
CONTINUE 

RETURN 

END 


SUBROUTINE TO CALCULATE CONSTANTS USED IN EVALUATION OF SOLID PHASE 
TEMPERATURES 


SUBROUTINE ONCNST (TSO, MA, CPBAR, UT, F BAR, KZC) 


COMMON 7 GEOTRY/R(45),EP (20) , DZ, THKCF (30) ,DV(20), OP, DPCF »M, 8, DT 
REAL Аа, ВНАС, МАСТ, PHC MACI I, RHCil MAE PHE , MAIO). CPZ CPCI 
BEZIECRECIT,CRE,UZEI 105,10) ,UET (10510) „UETITI (10, 10), UF (19716) 
FEAR RZCI]10, 10), BREI 010, 19), K8E (10,10), PI 

FE P CPBAP”(10) , UT (10, 10) EOIIE (19,109, ETWO (10,719), ETAP (10,10) 
EEAL EBAR (10, 10) , TSO (10, 10) 

REAL PR1,PR2,PR3, PR4 

PARAMETER {PI=3.14159,PB1=.000250,PP2=.0002,PP3=.00015) 
PARAMETER (PR4=.000117, PHZC=6300.,PHCI=1900., PHCII=1000.) 
PARAMETER (BHF210500.,CP2-2200.,CPCI-3000.,CPCII-3000.,CPFz200.) 
PARAMETEP. (KZ2C=40.,KC1=3.0,KCIl=1.5,KF=30.) 


DO 60 І-2,М-1 


CALCULATE MASS PER OMIT OUTSIDE SURFACE AREA 


MAZ=(PRIY*3-PP2**3) *BHZC/ (3* (PR1* %2)) 
MACI= (PA2%*3-PP34%3) *BHCI/ (3* (PP1**2)) 
MACII= (PR3®+3-PR4@# 3) "BHCII/ (3* (PAT**2)) 
MAF=(PR4**3) *PHF/ (3* (PR1**2)) 

MA (I) -МА2 + МАСТ+ МАСТТ + МАЕ 

ВЕ (г EO. 2) 45 (2) =. 0072 

EESTI EO. М-Т) MA(M-1)-z5.73 


CALCULATE AVEPAGE SPECIFIC HEAT 


CPBAR (I) =( (MAZ*CPZ+MACI *“CPCI+MACII*CPCII+MAF* CPF) /МА (Т)) 
IEZ(T 7.820222) CEBAR (2) —=422. 
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ТЕ (I ЕО М-1) СРВАВ (М-1) -155. 


DO 50 J=1,N 
KF (I,J) =1.7307*26.0* (AMAX1 (TSO (I, J) *1.8-459.67, 260.0) ) 
+ **(-0.1093) 
KZC(I,J)=22.67+.00867*TSO(I,J) 
КЕШГЕ, 7) АМАТ 7307 *(17921=19.7 *ALOG (1. 8*TSO (1. J) 


т -459.67)),10.0) 
IF (I .EQ. 2) THEN 
РТ. 


ЕВАВ (2,9) =1. 
KZC(2,J)=15. 
GOTO 50 
ENDIF 
IF (I .EQ. M-1) THEN 
UT(M-1,J)=1. 
FBAR(M-1,J)=1. 
KZC (M-1,J)=47. 
GOTO 50 
ENDIF 
“ALCULATE INTERNAL PARTICLE HEAT TRANSFER COEFFICIENTS 
UZC (I, J) =(PR2*PR1/ (PR1-PR2) ) *KZC (I,J) / (PR1**2) 
UCI (I, J) =(PR3*PR2/ (PR2-PR3)) *KCI (I,J) / (РВ1**2) 
UCII (I,J) =(PR4*PR3/ (PR3-PR4) ) *KCI (I,J) /(PR1**2) 
UF (I,J) =(2*KF (1I,J))/ (3*PR4) 
TOT ие A UCr (id) ) 4 (1 / UCI Lara) 
as О) 


CALCULATE THE SINGLE NODE FRACTIONS FOR TEMPERATURE DISTRIBUTION 
PONE (1 ,J)=0T (1,5) /UZC(1,J} 
ЕИО (Т, а) =ЕОМЕ (Т, с) + (0т (Т, с) UCI (TO) 
FTHR(L,J)=EFTWO(LI,J)+(UT(L,J)/UCILI(I;3)) 
FBAR(I,J)=(1/ (2*MA(I)*CPBAR(I)))*((MAZ*CPZ*FONE(I,J)) 
Б + (MACI*CPCI* (FONE (I,J) +FTWO(I,J))) 
+ + (MACII*CPCII* (FTHR(1,J)+FTWO(1,J)))+(MAF*CPF 
+ ETETHR(CI,J)GT1))) 
CONTINUE 
CONTINUE 
THE FOLLOWING MATERIAL PROPERTIES ARE ASSUMED FOR THE COLD FRIT 
DENSITY=7.95 E03 kg/m**3 
SPECIFIC HEAT=422 J/kg K 
CONDUCTIVITY=15 W/m K 
THE FOLLOWING MATERIAL PROPERTIES ARE ASSUMED FOR THE HOT FRIT 
DENSITY=2.10 E04 kg/m**3 
SPECIFIC HEAT=155 J/kg K 
CONDUCTIVITY=47 W/m K 
RETURN 
END 
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SUBROUTINE TO ADVANCE SOLID PHASE TEMPERATURES ONE TIME STEP 


SUBROUTINE SOLID (MA, CPBAR, UT, FBAR, TSO, TGO, PO, RHO, Q, VRO, VZO, 
+ KZC, TSST, UBAR, AV, ANN) 
COMMON/GEOTRY/R (45), EP (20), DZ, THKCF (30) , DV(20) , DP, DPCF, M, NH, DT 
COMMON/ BOUCON/PBI, PBO, TBI, RHOBI, HBI 
COMMON/ SOLVER/SOLSLD (100, 100) 
INTEGER I,J,K,M,N, UNITS 
REAL KG(10,10),KAP(10,10),PHI, PHII, PH, KEO, KRO, EMM, EDP 
REAL DP, DPCEF, KR, RE, THI, THII, CTHI,CTHII, MA(10) , CPBAR (10) 
REAL UT (10, 10) 
REAL FBAR(10,10),TSO(10,10),TSST(10,10), TGO (10, 10) , KSOL (10, 10) 
REAL AV(10),EP,VIS,GSO,FJH, TFLM, CPB, CPFLM, VISFLM, KFLM 
REAL-HLOG(0 10) PO(10,10),0(10.160), VBRO(X11, 10). V20(10, 11) 
REAL KW(10,10),KE(10,10),KS (10, 10) , KN (10, 10) 
REAL ACw( 0,10) ,ACE (10, 10), ACN(10/10) , AGS (10, 10) , ACO (10, 10} 
REAL CONVC (100) ,BC(10,10),UBAR(10,10) , KZ, RHO(10, 10) 
REAL ANN (10,10) ,KZ2C (10, 10) 
PARAMETER (PI-3.14159,UNITS-1,EMM-.78) 
INITIALIZE THE SOLUTION MATRIX 
DO 2 I=1,M*N 
DO 1 J=1,M*N 
SOLSLD (I, J) =0. 
CONTINUE 
CONTINUE 


CALCULATE GAS AND SOLID CONDUCTIVITIES 


THI=.68 
CTHI=.565 
THII=.144 
CTHII=.925 


DO 20 Т=2,М-1 
DO 10 J=1,N 
CALCULATE THE GAS PHASE CONDUCTIVITY 
KG €1 70) =PITCONDPO (1,90), 160 (1, J), UNITS) 
CALCULATE THE THERMAL CONDUCTIVITY IN THE SOLID 
KZ=KZC (I,J) 
IF (I .EO. M-1) THEN 


КЕО-К2 
GOTO 5 
ENDIF 


KAP (I, J) =KZ/KG (I,J) 
PHI=.5*((( (KAP (I,J) -1) /KAP (I,J) ) **2) *THI) / (ALOG (KAP (I, J) 


+ ~( (KAP (I, J)-1) *CTHI) ) - (KAP (I, J) -1) * (1-CTHI) /KAP (I, J)) 
1 =(27 3) “(17 (KAP (1, J5)))) 

PHIT=.5* ((( (KAP (I, J) -1) /KAP (I, J) )**2) *THII) / (ALOG (KAP (I, J) 
T -( (KAP (I, J)-1) *CTHII) ) -(KAP (I,J) ~1) * (1-CTHII) /KAP (I, J)) 


S27 3) "(1/7 ( KAP (1,5) ) 

PH=PHII+( (PHI-PHITI) * (EP (I) -.260) /.216) 

KEO=KG (1,3) * (EB (1) +( (1-EP (1) ) 7 (PH+ ( (273) *KAP (1,0) )))) 
IE (LEO. 2) THEN 
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EDP=DPCF ELSE 
EDP=DP 
ENDIF 
CALCULATE THE RADIATION CONTRIBUTION TO SOLID CONDUCTIVITY 
KRO=0..229*EMM*EP (I) *EDP* (TSO(I,J) **3)/ (10**6) 
KR=(1-EP (I)) /((1/KZ)+(1/KRO))+EP (I) *KRO 
COMBINE THERMAL AND RADIATION CONDUCTIVITIES TO GET OVERALL VALUE 
KSOL (I, J) =KEO+KR 
CALCULATE HEAT TRANSFER COEFFICIENT FOR SOLID TO GAS HEAT TRANSFER 
CALCULATE THE RE 
AV (I)23* (1-EP (I) ) / (EDP/2) 
IF (I .EQ. M-1) AV(M-1)=3659. 
WIS=PTVISCX(PO (lw), TGO( 1.) UNITS) 
GSO-SQRT (AMAX1 (ABS ( ( (VRO(I,J) *VRO(I*1,J))/2) **2) +ABS ( ( ( 
+ VZ0 (I,J) +VZ0 (I, J+1)) /2) **2), .00001)) *RHO (I,J) 
RE=GSO/ (AV(1)*VIS) 
DETERMINE RELATION FOR THE COULBURN FACTOR ON THE RE 
IF (RE .LT. 50.) THEN 
FJH=0.91* (RE** (-.51)) 
ELSE 
FJH=0.61* (RE** (-.41)) 
ENDIF 
TFLM=.5* (TSO(1,J)+TGO(1,J)) 
CPB=PTCP (PO(I, J), TGO(I,J), UNITS) 
CPFLM=PTCP (PO(I, J), TFLM, UNITS) 
VISFLM=PTVISC (PO(I,J), TFLM, UNITS) 
KFLM=PTCOND (PO (I,J), TELM, UNITS) 
HLOC (I,J) =FJH*CPB*GSO/ (( (CPFLM*VISFLM) /KFLM) ** (2/3) ) 
CALCULATE UBAR 
UBAR FOR THE FRITS IS ASSUMED TO BE EQUAL TO HLOC 
IF (I .EQ. 2) THEN 
UBAR (2, J) =HLOC (2, J) 
GOTO 9 
ENDIF 
IF (I .EQ. M-1) THEN 
UBAR (M-1, J) =HLOC (M-1, J) 
GOTO 9 
ENDIF 
DBAR(ISJ)EUTÉIÍ,J)*HLOC(IjJd)Z(EBAR(I,JG)*HLOC(I, 3) +UT (1, d))} 
IF (I „ЕО: 2) UBAR(2 J)=HLOC(I,J) 
CONTINUE 
WRITE (3, *)1,KSOL(I, J) , KRO, KR, KEO, HLOC (I,J), UBAR(I,J), 
+ AV(I),RE, PHI, PHII, PH, GSO 
CONTINUE 
CONTINUE 
CALCULATE EAST AND WEST FACE CONDUCTORS AND COEFFICIENTS BASED ON 
THE HARMONIC MEAN OF THE CONDUCTIVITIES ON EITHER SIDE OF THE 
INTERFACE 
DO 40 I=2,M-1 
DO 30 J=2,N 
KW(I,J)=2*KSOL (I,J) *KSOL (I, J-1) / (KSOL(I, J) +KSOL (I, J-1) ) 
KE (I1, J-1) =KW (I, J) 
ACW(I,J)-KW(I,J)* (R(2*I-1) **2-R(2* I*1) * *2) *PI/DZ 


169 











30 


40 
С 
C 


50 


60 


С 


С 


70 
80 
С 


90 
100 
C 


C 


ACE (I, J-1) =ACW(I, J) 
CONTINUE 
KW(I,1)=0. 
ACW(I,1)=0. 
KE(I,N)=0. 
АСЕ(Т Му-0. 
CONTINUE 
CALCULATE NORTH AND SOUTH FACE CONDUCTANCES AND COEFFICIENTS 
BASED ON THE HARMONIC MEAN CONDUCTIVIT 
DO 60 J=1,N 
DO 50 I=3,M-1 
F=(R(2*I-2)-R(2*I-1))/(R(2*I-2)-R(2*I)) 
KS(I,J)-(((1-F)/KSOL(I,J) ) F(F/KSOL(I-1,J))) ** (-1) 
KN (I-1, J) =KS (I, J) 
ACS (I,J) =KS (I,J) *2*PI*R(2*I-1) *DZ/ (R(2* 1-2) -R(2*I)) 
ACN (I-1,J)=ACS (I,J) 
CONTINUE 
KS(2,J)=0. 
ACS (2,J)=0. 
KN (M-1, J) =0. 
ACN (M-1,J)=0. 
CONTINUE 


CALCULATE THE CENTRAL POINT COEFFICIENT 
DO 80 I=2,M-1 
DO 70 J=1,N 
ACO (1, J) == (ACW(1, J) 4ACE (1,0) АСЗ СГ, и) +ACN (1,3) 
CALCULATE TIME DERIVATIVE COEFFICIENT 
BC (I, J) =MA (I) *CPBAR (I) *DV (I) *AV (I) / (DT*EP {I)) 
ANN (I, J) =(BC(I, J) -ACO (I,J) +UBAR (I,J) * (DV(I) /EP (I) ) *AV(I)) 
CONTINUE 
CONTINUE 
BUILD THE SOLUTION MATRIX AND CONSTANT VECTOR 
DO 100 I=2,M-1 
DO 90 J=1,N 
K= (1-2) *N+J 
SOLSLD (K, K) =-ANN (I,J) 
ТЕ (1 ІТ. М) SOLSLD(K) K41)=ACE (1, 2) 
IF (K*N .LE. (M-2)*N) SOLSLD(K, K+N) =ACN (I,J) 
IF (J .GT. 1) SOLSLD(K, K-1) =ACW(I,J) 
IF (K-N .GT. 0) SOLSLD(K, K-N) =ACS (I,J) 
CONVC (K) =-Q (I,J) -UBAR (I,J) * (DV( I) /EP(I) ) *AV(I) *TGO(I,J) 
+ -BC (I,J) *TSO(I,J) 
CONTINUE 
CONTINUE 
CALL SOLUTION ROUTINES 
NP=100 
CALL LUDCMP (SOLSLD, ((M-2) *N),NP, INDX, D) 
CALL LUBKSB(SOLSLD, ((M-2) *N),NP, INDX, CONVC) 
UPDATE THE SOLID PHASE TEMPERATURES 
I-2 
J=1 
DO 110 K=1, ( (M-2) *N) 
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TSST(1,J)=CONVC (K) 
IF ((K-(I-2)*N) .EQ. N) THEN 
J=1 
I=I+1 
ELSE 
J=J+1 
ENDIF 
CONTINUE 
TMCON1=MA (3) *CPBAR (3) /UBAR (3,1) 
TMCON2=MA (3) *CPBAR (3) /UBAR (M-2, 3) 
WRITE (6, *) ’ TMCON1’ , TMCON1, ' TMCON2' , TMCON2 
RETURN 
END 


SUBROUTINE TO PERFORM THE AXIAL MOMENTUM BALANCE 


SUBROUTINE AXLMOM (VR, VZ,VZO, RHO, RHOO, P, TG, BZAF, VZST,BVZ,AZF, 
+ AWZ,AEZ,ASZ,ANZ) 
COMMON/GEOTRY/R (45), EP (20), DZ, THKCF (30) ,DV(20),DP,DPCF,M,N, DT 
COMMON / BOUCON/PBI, PBO, TBI, RHOBI, HBI 
COMMON/ SOLVER/SOLZ (100, 100) 
INTEGER M,N,I,J,K,UNITS 
Reo VR(11, 10), V2(10,11),R, EP, RHO(10,10),P 110.10), TG (10/16) 
REAL MEZ(10,11),MW2(10,11),*582(10,11),MSZ (10,11), V28T (10, 11) 
EESLOAEZ(10 11),AW2(10,11),ANZ(10,11),ASZ(10, 11), AOZ 110, 11) 
ВЕАЬ В2 (10,11), АГЕ, АГЕВ, АГЕС, ALFD, GAM, GAMB, GAMC, $Z (10,11) 
REAL DT,DZ, VSR, THKCF, DVZ, SOLZ, CONVZ (100), BZAF (10,11), ARI, ARO 
REAL FRIC(10,11),REI,REO, VISI, VISO, HDI, HDO, BVZ(10,11),AZF (10,11) 
REAL VZ0(10,11),RHOO(10,10),BZ0(10, 11), RHOOS 
PARAMETER (PI=3.14159, UNITS=1) 
INITIALIZE THE COEFFICIENT MATRIX 
DO 8 I=1,M*N 
DO 4 J=1,N*M 
SOLZ (I, J)=0. 
CONTINUE 
CONTINUE 
GENERAL FORM OF THE EQUATION: 
-(BZ-AOZ)*VZO 4 AEZ*VZE + AWZ*VZW + ANZ*VNZ + ASZ*VZS = -SZ - BZ*VR 


CALCULATE CONTROL VOLUME FACE MASS FLUXES AND VELOCITY COEFFICIENTS 


SOUTH AND NORTH 
DO 16 J=2,N 
DO 12 I=2,M 
UPWIND FACTORS FOR DENSITY 
IF (VR(I,J-1) .GE. 0.) THEN 
ALFB=0.5 
ELSE 
ALFB=-0.5 
ENDIF 
тр Вс) GE) 0.) THEN 
ALFC=0.5 
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ЕТСЕ 
ALFC=-0.5 
ENDIF 
C FACE FLUX 
MSZ (I,J) =(((0.5+ALFB) *RHO (I-1,J3-1)+(0.5-ALFB) *RHO(I,J-1)) 
*VR(I,J-1)+((0.5+ALFC) *RHO(I-1,J)+(0.5-ALFC) *RHO(I,J)) 
n *VR(I,J))*PI*R(2*I-1)*(DZ) 
MNZ(I-1,J)-2 MSZ(I,J) 
C | UPWIND FACTOR FOR MOMENTUM 
IF ((VR(I,J-1)+VR(I,J)) .GE. 0.) THEN 
ALFD=0.5 
ELSE 
ALFD=-0.5 
ENDIF 
C VELOCITY COEFFICIENT 
ASZ (I,J) =MSZ (I,J) *(0.5+ALFD) 
ANZ(I-1,J)-2-MNZ(I-1,J)*(0.5-ALFD) 
на CONTINUE 
C | BOUNDARY CONDITIONS 
MSZ(1,J)-0. 
non 29. 
ММ2 (М, Ј) =0. 
АМ2 (М, Ј) =0. 
16 СОМТІМОЕ 
MSZ (1,1) =0 
ASZ (1,1) =0 
TE (VR(2, l) .GE. 0.) THEN 
ALFC=0.5 
ELSE 
ALFC=-0.5 
ENDIF 
MNZ (1,1)=((0.5+ALFC) *RHO(1,1)+(0.5-ALFC) *RHO(2,1))*VR(2,1)*PI 
+ * (DZ) *R(3) 
АМ№2 (1,1) =-ММ7 (1,1) * (0.5-ALFC) 
ІЕ (ҮБ(М,1) „СЕ. 0.) THEN 
ALFC=0.5 
ELSE 
ALFC=-0.5 
ENDIF 
MSZ (M, 1)2((0. 5*ALFC) *RHO(M-1, 1) * (0. 5S-ALFC) *RHO (M, 1) ) *VR (M, 1) 
+ *PI*R(2*M-1) * (DZ) 
C EAST AND WEST FACE FLUXES AND COEFFICIENTS 
DO 28 I=1,M 
DO 24 J=2,N 
MWZ (I,J)zRHO(I,J-1)*((VZ(I,J-1)4VZ(I,J))/2) * (R(2*I-1) **2 
+ -R(2*I+1) **2) *PI 
MEZ (I,J-1)=MWZ (I, J) 
C | UPWIND FACTOR FOR MOMENTUM 
ЕССЕ ОУ ТУСТЕ су) 4 GE-20-) ^ THEN 
GAMB-0.5 
ELSE 
GAMB--0.5 
ENDIF 
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VELOCITY COEFFICIENT 
AWZ (I1, J) =MWZ (I, J) * (0. 5+GAMB) 
AEZ (I, J-1) =-MEZ (I, J-1) * (0.5-GAMB) 
CONTINUE 
CONTINUE 
BOUNDARY CONDITIONS 
DO 29 1=2,M-1 
IF (VR(I,1) .GE. 0.) THEN 


ALFC=0.5 
ELSE 
ALFC=-0.5 
ENDIF 
MSZ2(T,1)=((0.3+ALEC) *RHO(I=1,1)+(0.5-ALFC)*RHO(1,1)) 
+ PURTA s PI#B (2151) DE 


MNZ (I-1,1)=MSZ (1,1) 
ASZ (1,1)=MSZ(1,1)*(0.5+ALFC) 
ANZ (I-1,1)=-MNZ (I-1,1)*(0.5-ALFC) 
MWZ (I,1)=0. 
AWZ(I,1)-0. 
CONTINUE 
BOUNDARY FLUXES AND COEFFICIENTS 
EAST FACES NOT CO-LOCATED WITH WEST FACES 
DO 32 I-1,M 
MEZ (I,N)-RHO(I,N)*(VZ(I,N)/2) *PI* (R(2*I-1) **2-R(2*I41) 
+ жж2) 
IE (VEI NY) „СЕ. 0.) THEM 
АЕ? (Т,М) =0. 
ELSE 
AEZ (I,N)--MEZ(I,N) 
ENDIF 
CONTINUE 
INLET AND OUTLET TO PLENUMS 
MWZ (1,1) =RHOBI*VZ (1,1) * (R (1) **2-R (3) **2) *PI 
MWZ (M, 1) -RHO (M, 1) *VZ (M, 1) * (R(2*M-1) **2) *PI 
AWZ (1,1)=0. 
AWZ (M, 1) =0. 
CALCULATE NODE CENTRAL POINT COEFFICIENT 
DO 44 I=1,M 
DO 40 J=1,N 


AOZ(I,J)=-(AEZ(I,J)+AWZ(I,J)+ANZ(I,J)+ASZ(I,J)+MEZ(I,J) 
+ -MWZ (I, J) +MNZ (I,J) =MSZ2 (I, J} ) 
CONTINUE 
CONTINUE 


CALCULATE CONSTANT TERMS 
DO 52 J=2,N 
DO 48 I=1,M 
UPWIND FACTOR FOR DENSITY 
ТЕ (VZ(I,J) .GE. 0.) THEN 
GAM=0 .5 
ELSE 
GAM=-0.5 
ENDIF 
RHOS=(0.5+GAM) *RHO (I, J-1)+(0.5-GAM) *RHO (I, J) 
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RHOOS=(0.5+GAM) *RHOO (1, J-1)+(0.5-GAM) *RHOO (1, J) 
DVZ=PI* (R(2*I-1) **2-R(2*I+1) **2) *DZ 
BZ(I, J) =RHOS*DVZ*EP (I) /DT 
BZO(1,J)=RHOOS*DVZ*EP (1) /DT 
BVZ (I, J)=DVZ*EP (I) /DT 
SUPERFICIAL VELOCITY 
VZS=(((.25* (VR(I+1, J) +VR(I+1,J-1) +UR (I, J}4VRUL, T-1)))**2 
+ +VZ (I,J) **2)**0.5) 
CORRECTION FOR COLD FRIT VARIABLE THICKNESS 
IF (I .EQ. 2) DVZ=(THKCF (J) * THKCF (J-1) ) *PI*R(2*I)*DZ 
SOURCE TERM 
SORC- (P (I,J-1) -P(I,J)) * (DVZ/DZ) *EP (I) 
FRICTION TERM 
PS-(0.54GAM) *P(I,J-1) 4 (0. 5- GAM) *P (I,J) 
TS- (0.54GAM) * TG(I,J-1) * (0.5-GAM) * TG (I, J) 
VIS=PTVISC (PS, TS, UNITS) 
IF (I .EQ. 2 .OR. I .EQ. (M-1)) THEN 
DPA-DPCF 
ELSE 
DPA-DP 
ENDIF 
FRIC(I,J)-DVZ*((150*VIS*((1-EP(I)) **2)/ ((DPA**2) * 
E (EP (I) **3)))+(1.75*RHOS*ABS (VZS) *(1-EP (I))/ (DPA* 
T (ЕРТ КЕЗ EPI) 
SZ (I,J) =SORC 
WRITE (3, *)I,J,DPA, EP (I), DVZ 
WRITE(3,*)SORC, V2S, VZSUP, VIS 
CONTINUE 
CONTINUE 
INLET AND OUTLET PLENUMS 
SOURCE TERMS 
ARI=PI* (R(1)**2-R(3)**2) 
ARO=PI* (R(2*M-1)**2) 
SZ(1,1)=(PBI-P(1,1))*ARI 
SZ (M, 1) - (PBO-P (M, 1) ) *ARO 
BZ(1,1)-RHOBI*ARI*DZ/DT 
BZO(1,1)-BZ(1,1) 
BZ (M, 1) -RHO (M, 1) *ARO*DZ/DT 
BZO (M, 1) xRHOO (M, 1) *ARO*DZ/DT 
BVZ (1,1)=BZ(1,1)/RHOBI 
BVZ (M, 1) =BZ (M, 1) /RHO (M, 1) 
FRICTIONAL LOSSES IN THE PLENUMS 
CALCULATE HDI, HDO, THE HYDRAULIC DIAMETERS 
HDI=2*(R(1)-R(3)) 
HDO=2*R (2*M-1) 
CALCULATE REYNOLDS NUMBERS, UPWINDING OF DENSITIES CONSIDERED 
UNNECESSARY 
DO 53 J=1,N 
VITIST=PTVISC(P (1, J), Te (1, J), UNITS) 
REI=RHO(1,J)*ABS(VZ(1,J))*HDI/VISI 
VISO=PTVISC(P (M, J) , TG (M, J) , UNITS) 
REO-RHO (M, J) *ABS (VZ (M, J) ) *HDO/VISO 
CALCULATE FRICTION FORCE USING ACURVE FITTED TO THE MOODY DIAGRAM 
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FOR FRICTION FACTOR CALCULATIONS 
IF (REI .LT. 1000.) THEN 
FRIC(1,J)=.0482*RHO(1,J) *DZ*ARI/ (2*HDI) 
GOTO 531 
ENDIF 
FRIC(1,J)=.138* (REI**(-.151)) *RHO(1,J) *DZ*ARI*ABS (VZ(1,J}) 
+ / (2*HDI) 
IF (REO .LT. 1000.) THEN 
FRIC(M, J) -.0482*RHO (M, J) *DZ*ARO/ (2*HDO) 
GOTO 53 
ENDIF 
FRIC(M, J)=.138* (REO** (-.151)) *RHO (M, J) *DZ*ARO*ABS (VZ (M, J)) 
+ / (2*HDO) 
CONTINUE 
DO 55 I=1,M 
BO. 54 J=2,N 
BZAF (I, J) =BZ (I,J) -AOZ(I,J)+FRIC(I,J) 
AZF (I,J) =- (BZAF (I,J) -BZ(I,J)) 
WRITE (3,*) BZ(I,J),FRIC(I,J), BZAF (I,J) 
CONTINUE 
CONTINUE 
BZAF (1,1) =BZ (1,1) -AOZ (1, 1) +FRIC (1,1) 
BZAF (M, 1) =BZ (M, 1) -AOZ (M, 1) *FRIC (M, 1) 
AZF (1,1) =- (BZAF (1,1) -BZ (1,1) ) 
AZF (M, 1) =- (BZAF (M, 1) -BZ (M, 1) ) 
WRITE(3- *)BZAE(L,1),BZ(1,1),SZ(1,1), FRIC (T, 1) 
WRITE(3,*)BZAF (M,1),BZ (M, 1), SZ (M, 1) , FRIC (M, 1) 
FORCE WALL VZ TO ZERO 
DO 56 I=2,M-1 
SZUT,1)-0. 
RHOS=RHO (I, 1) 
RHOOS=RHOO (I, 1) 
DVZ=PI* (R(2*I-1) **2-R(2*I+1) **2) *DZ 
BZ (I, 1)=RHOS* (DVZ/2) *EP (I) /DT 
BZO (I, 1)=RHOOS* (DVZ/2) *EP (I) /DT 
BVZ (I,1)=DVZ*EP (I)/ (2*DT) 
VZS=.5*(VR(I,1)+VR(I+1,1)) 


DPA=DPCF 
PRIC (1,1) =bva*((150*ViS* ((l-Be( 2) ) ТОИ СОВА) А 
+ (BP (Lj <3) ) y+ (1. 7S*RHOS*ABS(VZS) *(1-EP(1))/ (DPA* 
Т (ЕР (1) **3)))) *ЕР(І) 


BZAr (172) =(BZ(1,1) -ACZ (1; 1)+FRIC(I, 1))) 
AZF(I,1)-AOZ(I,1)-FRIC(I,1) 
CONTINUE 
CONSTRUCT THE SOLUTION MATRIX AND CONSTANT VECTOR 
ENTERING THE MANIFOLD TURNING COEFFICIENT FOR THE INLET AND 
OUTLET PLENUMS 
CTD=1.05 
СТС=2.66 
DO 64 J=1,N 
SOLZ(J,Jy=- (BZ (1, J) -(CTD*AO2Z (1, J) ) +FRIC(1, J) ) 
ЕТ. м) SOBZ(J,UTI)CID*ABEZ(L;«) 
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LE (J GI. 1) SOLA (J,.0-1) =CTh*Awz (1,7) 
SOLZ (J, J+N) =ANZ (1, J) 
CONVZ (J) 2-SZ2(1,J) -(BZO(1, J) *VZO(1,J)) 
DO 60 I-2,M-1 
ВЕ МЕЛ 
SOLZ (K, K) =-BZAF (I,J) 
ТЕ (0 .LT. N) SObLZ (K,K+1) =AEZ (1,0) 
IF (K+N .LE. ((M-2)*N)) SOLZ(K, K+N) =ANZ (I,J) 
IF (J ст. 1) 5002 (К,К-1) +АЙ? (ГТ, ч) 
ТЕ (К-М „СТ. 0) 5012 (К, К-Н) -АЗ2 (Г, 7) 
CONVZ(K)z-SZ(I,J)-(BZO(I,J)*VZO(I,J)) 
CONTINUE 
CONTINUE 


SKIP OVER HOT FRIT SINCE THERE IS NO AXIAL MOTION AND THIS WOULD 


RESULT IN A SOLUTION MATRIX WITH A ROW OF ZEROS WHICH IS 
INCOMPATIBLE WITH THE SOLUTION METHOD USED. 
DO 65 J-1,N 
K=(M-1) *N+J 
SOLZ (K, K) =- (BZ (M, J) - (CTC*AOZ (M, J) ) *FRIC (M, J) ) 
IF (J .LT. N) SOLZ(K,K+1)=AEZ (M, J) *CTC 
ЕТ СТУ) Зое (К КЕ) =AWz (Mpg) *CTC 
CONVZ (K) =-SZ (M, J) - (BZO (M, J) *VZO (M, J) ) 
CONTINUE 
NP=100 
DO 66 K=1,45 
WRITE (3, *) K, SOLZ(K, K) 
CONTINUE 
CALL SOLUTION SUBROUTINES 


NOTE THET THE UPDATED VELOCITIES ARE RETURNED AS A MODIFIED CONVZ 


VECTOR 
CALL LUDCMP (SOLZ, (M*N),NP, INDX,D) 
CALL LUBKSB (SOLZ, (M*N),NP, INDX, CONVZ) 
CONVERT CONVZ TO I,J FORM 
I=1 
J=1 
DO 68 K=1, (M*N) 
VZST (I, J) =CONVZ (K) 
IF ((K= (I=1) *N) -ЕО. МЫ) THEN 
J=1 
I=I+1 
ELSE 
J=J+1 
ENDIF 
CONTINUE 


* 


CORRECTIONS TO CONSTANTS USED IN LATER PARTS OF THE PROGRAM TO 


REFLECT THE TURNING COEFFICIENTS 
DO 69 J=1,N 

AZF (1, J) =A0Z (1, J) *CTD-FRIC (1, J) 

AZF (M, J) =AOZ (M, J) *CTC-FRIC(M, J) 

B2ZAF (1, J) =(BZ (1,J) - (AOZ (1, J) *CTD) +FRIC (1, J) ) 
BZAF (M, J) = (BZ (M, J) - (AOZ (M, J) *CTC) +FRIC (M, J) ) 
AEZ (1,J) -AEZ (1, J) *CTD 

AWZ (1,J)2AWZ(1,J) *CTD 
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AEZ (M, J)=AEZ (M, J) «СТС 
AWZ (M, J) =AWZ (M, J) *CTC 
CONTINUE 
RETURN 
END 


SUBROUTINE TO PERFORM THE RADIAL MOMENTUM BALANCE 
MODIFIED TO WORK IN SUPERFICIAL VELOCITIES 


SUBROUTINE RADMOM (VR, VRO, VZ, RHO, RHOO, P, TG, BRAF, VRST, BVR, 

+ ARF, AWR, AER, ASR, ANR) 
COMMON/GEOTRY/R (45), EP (20) , DZ, THKCF (30) ,DV(20), DP, DPCF,M,N, DT 
COMMON/BOUCON/PBI,PBO, TBI, RHOBI, HBI 
COMMON/SOLVER/SOLR (100, 100) 
INTEGER M,N,I,J,K,UNITS 
REAL VR(11,10),V2(10,11),R,EP,REO(10,10),P (10,10), TG(LO, 10) 
REAL MER(11,10),MWR(11,10),MNR(11,10),MSR(11,10), VRST(11, 10) 
REAL AER(11,10),AWR(11,10),ANR(11,10),ASR(11,10) 
REAL AOR(11,10),BR(11,10),ALF,ALFB, ALFC, ALFD, GAM, GAMB, GAMC 
REAL ЗВ (11,10), ЕВТС, ВВАЕ (11,10), ВУВ (11,10), АВР (11,10) 
REAL DT, DZ, VSR, THKCF, DV, SOLR, CONVR (100) 
КЕАІ УВО (11,10), ВНОО (10,10) ,ВВО (11,10), ВНООЅ 
PARAMETER (PI=3.14159, UNITS=1) 


GENERAL FORM OF THE EQUATION: 
— (BR-AOR) *VRO + AER*VRE + AWR*VRW + ANR*VRN + ASR*VRS = -S - BR*VR 


CALCULATE CONTROL VOLUME FACE MASS FLUXES AND VELOCITY COEFFICIENTS 


INITIALIZE SOLR 
DO 2 I=1,M 
DO 1 J=1,N 
SOLR(I, J) =0 
CONTINUE 
` CONTINUE 
SOUTH AND NORTH 
DO 10 I=2,M 
DO 8 J=1,N 
SOUTH 
ESTIMATE VELOCITY AT THE INTERFACE USING: 
VRi*AREAi-VR* (R(2I-1)/R(21-2))*R(21-2) *2*PI*DZ 
MSR(I,J)-RHO(I-1,J) *VR(I, J) *R(2*I-1) *2*PI* 
+ DZ 
NORTH 
MNR (I-1,J)=MSR(I,J) 
IF ((VR(I,J)*VR(I-1,J)) .GE. 0.) THEN 
ALFB=0.5 
ELSE 
ALFB=-0.5 
ENDIF 


C COEFFICIENTS 
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20 


ASR(I,J)=MSR(I,J)*(.5+ALFB) 
ANR (I-1,J)=-MNR (I-1,J)*(.5-ALFB) 
CONTINUE 
CONTINUE 
BOUNDARIES 
DO 12 J=1,N 
MNR (M, J) =RHO (M, J) *VR (M, J) *2*PI*R(2*M) *DZ 
ANR (M, J)=0.0 
CONTINUE 
EAST AND WEST 
DO 20 I=2,M 
DO 18 J=2,N 
FLUX UPWIND FACTORS 
ТЕ оу СЕ 0.) ТИЕҢ 
GAMB=0.5 
ELSE 
GAMB=-0.5 
ENDIF 
ШЕ има т-1 т) СЕ. 0.) THEN 
GAMC=0.5 
ELSE 
GAMC--0.5 
ENDIF 
WEST FACE FLUX 
MWR(I,J)2((.54GAMB) *RHO(I,J-1)*(.5-GAMB) *RHO(I, J) ) *VZ(I,J) 


+ *PI* (R(2*I-1) **2-R(2*I)**2)+((.5+GAMC) *RHO (I-1,J-1)+ 
+ (0.5-GAMC) *RHO (I-1,J) ) *VZ{I-1, J) *PI* (R(2*I-2) **2-R (2*I-2) 
+ **2) 


EAST FACE FLUX 
MER (I, J-1) =MWR (I,J) 
MOMENTUM UPWIND 
ТЕО СТЕСТ О СЕ ОС ЕНЕМ 
GAMBC-0.5 
ELSE 
GAMBC--0.5 
ENDIF 
WEST VELOCITY COEFFICIENT 
AWR (I, J) =MWR (I, J) * (0. 5+GAMBC) 
EAST VELOCITY COEFFICIENT 
AER (I, J-1) =-MER (I, J-1) * (0. 5-GAMBC) 
CONTINUE 
CONTINUE 
BOUNDARY CONTROL VOLUMES 
VZ (1,1) ASSUMED GE O AND VZ(M,1) LT O 
MWR (2, 1) =RHOBI* (R (2) **2-R (3) **2) *PI*VZ (1,1) 
AWR (2, 1) =MWR (2, 1) 
MWR (M, 1) =RHO(M, 1) *PI* (R(2*M-1) **2-R(2*M) **2) *VZ (M, 1) 
AWR (M, 1) =0. 
MER (2,N) =0. 
AER (2,N)=0. 
MER (M,N)=0. 
AER (M,N)=0. 
DO 24 I=3,M-1 
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28 
32 


MWR (I,1)=.0 


АНЕ (1,1) =0. 

МЕК (І, №) =0. 

AER(I,N)-20. 
CONTINUE 


CALCULATE THE CENTRAL POINT COEFFICIENT 
DO 332 ТМ 
DO 28 J=1,N 
AOR (I,J)=- (AER(I,J) +AWR(I,J)+ASR(I,J)+ANR(I,J)+MER(I, J) 
+ -MWR (I,J) +MNR(I,J)-MSR(I,J)) 
CONTINUE 
CONTINUE 
CALCULATE THE CONSTANT TERMS 
DO 40 J=1,N 
DO 36 I=2,M 
IF (VR(I,J) .GE. O.) TEEN 
SER 05 
ELSE 
ALF=-0.5 
ENDIF 
CALCULATE GAS VOLUMES USED TO CALCULATE BR(I,J) 
DVA- (R(2*I-2) **2-R(2*I-1)**2) *PI*DZ*EP (I-1) 
DVB= (R(2*I-1)**2-R(2*I)**2) *PI*DZ*EP (I) 
RHOS=( (0.5+ALF) *RHO(I-1,J)+(0.5-ALF) *RHO(I,J)) 
RHOOS= ( (0.5+ALF) *RHOO (I-1,J) +(0.5-ALF) *RHOO(I,J)) 
FIND SMALLER VOID FRACTION IN THE CONTROL VOLUME 
EPS=MIN (EP (I-1),EP(I}) 
BR(1,J)=RHOS* (DVA+DVB) /DT 
BRO (I,J) =RHOOS* (DVA+DVB) /DT 
BVR(I,J)=BR(I,J)/RHOS 
SUPERFICIAL VELOCITY WEIGHTED BY VOLUMES 
VRSB=((0.5*(VZ(I,J5) +VZ (I, J+1)) **2) +VRSC**2) **0.5 
VRSA=((0.5* (VZ (I-1,J) +VZ (I-1, J+1) ) **2) +VRSC**2) 
+ **0.5 
PS=(0.5+ALF) *P(I-1,J)+(0.5-ALF)*P(I,J) 
TIS=(0.5+ALF) *TG(I-1,J)+(0.5-ALF) *TG(I,J) 
VIS=PTVISC (PS, TIS, UNITS) 
CALCULATION OF THE SOURCE TERM 
AREAA=PI*D2* (R(2*I-2)+R(2*I)) 
PRESSURE DIFFERENCE ACTS ON SMALLEST GAS AREA OF THE CONTROL VOLUME 
SORC-AREAA* (P(I-1,J)-P(I,J)) *EPS 
FRICTION TERM 
THKA-R (2*I-2)-R(2*I-1) 
THKB-2R (2*I-1)-R(2*I) 
CORRECTION FOR THE VARIABLE THICKNESS OF THE COLD FRIT 
IF (I .EQ. 2) THEN 
THKB-THKCF (J) /2 
DPB-DPCF 
DPA-DP 
ELSE 
DPA=DP 
DPB=DP 
ENDIF 
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ТЕ (I .EQ. 3) THEN 
THKA=THKCF (J) /2 


DPA=DPCF 
DPB=DP 
ENDIF 
FRICA=THKA* ( (150*VIS* ( (1-EP (I-1) ) **2) / ( (DPA**2) * 
+ (EP (I-1) **3)))4(1.75* RHOS*VRSA* (1-EP (1-1) ) / (DPA* 
T (BET), 
FRICB=THKB* ( (150*VIS* ( (1-EP (I) ) **2)/ ((DPB**2) * 
+ (EP (I)**3))) * (1.75*RHOS*VRSB* (1-EP (I) ) / (DPB* 
t (EP (1) 4 *3)))) 


FRIC=(FRICA+FRICB) *AREAA*EPS 

SR (I,J) =SORC 

BRAF (I, J)=(BR(I, J) -AOR (I, J) +FRIC) 
ARF (I, J) =- (BRAF (I, J) -BR(I,J)) 


С WRITE (3, *) I, J, DPA, DPB, EP (І) 

€ WRITE (3, *) THKA, THKB 

С WRITE (3, *) SORC, FRICA, FRICB, FRIC 
C WRITE (3, *) VRSA, VRSB, VRSC, SR(I,J) 
36 CONTINUE 

40 CONTINUE 


C CONSTRUCT THE SOLUTION MATRIX AND CONSTANT VECTOR 
DO 48 I=2,M 
DO 44 J=1,N 
K= (1-2) *N+J 
SOLR (K, K) =-BRAF (I, J) 
IF (J .LT. N) SOLR(K, K+1) =AER (I,J) 
ТЕ ((K+N) .LE. (M-1)*N) SOLR(K, K+N) =ANR (I,J) 
IF (J СТО Т) SOLDR(KI K 1)=SAWR(I,J) 
ТЕ (K-N .GT. 0) SOLR(K, K-N) =ASR(I,J) 
CONVR (K) =-SR(I, J) -(BRO(I, J) *VRO(I,J)) 


44 CONTINUE 
48 CONTINUE 
NP=100 


C CALL SOLUTION SUBROUTINES 
NOTE THAT THE UPDATED VELOCITIES ARE RETURNED AS A MODIFIED CONVR 
Е VECTOR 
CALL LUDCMP (SOLR, ((M-1)*N), NP, INDX, D) 
CALL LUBKSB(SOLR, ((M-1) *N),NP, INDX, CONVR) 
a CONVERT UPDATED VELOCITIES TO I,J FORM 
I=2 
J=1 
DO S2 K=1 (M=1) *N 
VRST (I, J) =CONVR (K) 
ГЕ ((Е-1-2) “М) EQ. N) THEN 
J=1 
I=I+1 
ELSE 
J=J+1 
ENDIF 
92 CONTINUE 
RETURN 
END 


о 
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THE FIRST MOMENTUM CORRECTOR - THE PRESSURE EQUATION 


SUBROUTINE PRESUR (VR, VZ,RHO, RHOO, RHST,P,PO, PST, PIO, TG, 

+ BZAF, BRAF, VRSTT, VZSTT) 
COMMON/GEOTRY/R(45),EP(20),DZ, THKCF (30) ,DV(20) , DP, DPCE, M, N, DT 
COMMON/BOUCON/PBI, PBO, TBI, RHOBI,HBI 
COMMON/ SOLVER/SOLP (100, 100) 

GENERAL FORM OF THE EQUATION 
NOTE THAT THIS SUBROUTINE SOLVES FOR THE PRESSURE INCREMENT 
- (RHO*DV/ (DT*P)+DO)*PIO + DE*PIE + DW*PIW + DS*PIS + DN*PIN =DIVFLUX 


INTEGER M,N,1,J,K,UNITS, IPVT(100) 
PEAL УВ (11,10), \/2 (10, 11), ВНО (10,10), ВНЗТ(10,10),Р (10,10), 
WerG (10, 10), 302 (10,11), B2(10, 11), AOR (11,10),,B8 (11, 10) ,SOLP, 
+ CONVP (100), VRSTT(11,10),VZSTT(10,11),BRAF (11,10), BZAF (10, 11) 
REAL R,EP,DZ,DV 
REAL GAM, PST(10,10) 
READ DO (10,10), DE(LO, 10) DW(10), 10), Ds (10,10) DN(10, 10), 
+ ELUXE (10,10), FLUXW(10,10),FLUXS (10,10) ,FLUXN(10,10),PIO(10,10), 
+ DIVFLX(10,10) 
REAL RHOO (10,10),PO(10,10) 
PARAMETER (PI=3.14159, UNITS=1) 
INITIALIZE THE COEFFICIENT MATRIX 
DO 4 I=1,M*N 
DO 2 J=1,M*N 
SOLP (I,J) =0. 
CONTINUE 
CONTINUE 
CALCULATION OF THE FLUXES AND PRESSURE INCREMENT COEFFICIENTS 


EAST AND WEST FACES 
WEST BOUNDARY 
DO 8 I=2,M-1 
DW(I,1)=0. 
FLUXW(1,1)=0. 
CONTINUE 
INLET, ASSUMING THE VELOCITY IS POSITIVE 
DW(1,1)-RHOBI*(((R(1)**2-R(3) **2) *PI) **2) /BZAF (1, 1) 
FLUXW(1,1)=RHOBI*VZ (1,1) *PI*(R(1) **2-R(3) **2) 
OUTLET, ASSUMING VELOCITY IS NEGATIVE 
DW (М, 1) =RHO (M, 1) * (((R(2*M-1)**2)*P1)**2)/BZAF (M, 1) 
FLUXW (M, 1) -RHO (M, 1) *VZ (M, 1) *PI* (R(2*M-1)**2) 
EAST BOUNDARY 
DO 16 I=1,M 
DE(I,N)=0. 
FLUXE (1,N)=0. 
CO-LOCATED EAST AND WEST FACES 
рО 12 3=2,N 
UPWIND FACTOR FOR DENSITY 
IF мито GE. 0.) THEN 
GAM=0 .5 
ELSE 
GAM=-0.5 
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ENDIF 
WEST - EAST PRESSURE INCREMENT COEFFICIENTS 
DW(I, J) =((0.5+GAM) *RHO(I, J-1) +(0.5-GAM) *RHO (I,J) ) 
+ *((PI* (R(2*I-1) **2-R(2*I+1) **2)) **2) *EP (I) /BZAF (I,J) 
DE (I, J-1)=DW (I, J) 
FLUXW (I,J)=((0.5+GAM) *RHO (I, J-1) + (0.5-GAM) *RHO (I, J) ) 


+ Жуа (ІІ ТЕТТК(2%1-1) <%2-БЕ (27141) %%2) 
ELUXE (I, J-1) =FLUOXW (I, J) 
CONTINUE 
CONTINUE 


SOUTH AND NORTH FACES 
BOUNDARY CONDITIONS 
DO 28 J=1,N 
DN(M,J)=0. 
FLUXN (M, J) 20.0 
0$ (1,9) =0. 
FLUXS (1,J)=0.0 
CO-LOCATED NORTH AND SOUTH FACES 
DO 24 I=2,M 
UPWIND DENSITY 
ТЕ ((УЕ(І,7)) .GE. 0.) THEN 
ALF=0.5 
ELSE 


ENDIF 
EPS=MIN (EP (1) EP (1-1)) 
NORTH SOUTH PRESSURE INCREMENT COEFFICIENTS 
DS(1,J)=((0.5+ALF) *RHO(I-1,J)+(0.5-ALF) *RHO(1,J)) 
+ * ((2*PI*R(2*I-1)*DZ)**2)*EPS/(BRAF (1,J)) 
DN (I-1, J) =DS (I,J) 
FLUXS (I,J) =((0.5+ALF) *RHO (I-1,J)+(0.5-ALF)*RHO(I,J)) 


+ “VB (1, T) *2*PTI*R12"%*1=1)%D2 
FLUXN (I-1, J) =FLUXS (I,J) 
CONTINUE 
CONTINUE 


CENTRAL POINT COEFFICIENTS AND FLUX DIVERGENCE 
DO 36 I=1,M 
DO 32 J=1,N 
DC (I, J)=DE (I, J) +DW(I, J) +DS (I, J) +DN (I,J) 
DIVFLX (I, J) =FLUXE (I, J) -FLUXW (I, J) +FLUXN (I, J) -FLUXS (I, J) 
CONSTRUCT SOLUTION MATRIX AND CONSTANT VECTOR 
K=(I-1)*N+J 
SOLP (K, K) =~ ((DV(1) /DT) * (RHOO (I,J) /PO(I,J))+DC(I,J)) 
IF (J .LT. N) SOLP(K,K*1)-DE(I,J) 
ТЕ ((K+N) .LE. (M*N)) SOLP (K,K*N) -DN (I, J) 
IF (J .GT. 1) SOLP (K,K-1) =DW (1,3) 
IF ((K-N) .GT. 0) SOLP (K,K-N)=DS (1,3) 
CONVP (K) =DIVFLX (I, J) 
WRITE (3, *) 1, FLUXS (I,J) , FLUXN (I,J) , FLUXW(I, J) , FLUXE (I,J) 
WRITE (3 2)DS(”(L J) DNC(1 J) DW(T J) DB I1 3) DO(L 2) 
WRITE (3, *) DV(I) , RHOO(I,J),PO(I,J),BZAF (I,J), BRAF (I,J) 
CONTINUE 
CONTINUE 
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NP=100 
MP=1 
MPP=1 
NR=M*N 
CALL SOLUTION SUBROUTINES 
CALL LUDCMP (SOLP,NR,NP, INDX,D) 
CALL LUBKSB(SOLP,NR,NP, INDX, CONVP) 
SOLUTION RETURNED AS MODIFIED CONVP VECTOR 
UPDATE PRESSURES 
I=1 
J=1 
DO 40 K=1,M*N 
PST(I,J)=PO(I,J)+CONVP (K) 
PIO(I,J)=CONVP (K) 
WRITE (3,*)P(I,J),PIO(I,J) 
IF ((K-(I-1)*N) .EQ. N) THEN 
J=1 
Т=1+1 
ELSE 
J=J+1 
ENDIF 
CONTINUE 
UPDATE DENSITIES BASED ON NEW PRESSURES 
DO 52 1=1,M 
DO 44 J=1,N 
RHST (1,J)=PTDENS (PST(1,J),TG(I,J),UNITS) 
CONTINUE 
UPDATE VELOCITIES 
AXIAL VELOCITIES 
GENERAL FORM OF THE EQUATIONS 
VZSTT(I, J)=(DW(I,J) * (PIO(I, J-1) -PIO(I, J) )+VZ(1I,J) *RHZOLD*AREA) 
/ (RHZ* AREA) 
IF (I .EQ. M-1) GOTO 52 
DO 48 J=2,N 
UPWIND TO COMPUTE DENSITY AT VZ 
IF (VZ(I,J) .GT. 0.) THEN 
RHZOLD-RHO(I,J-1) 
RHZ-RHST(I,J-1) 
ELSE 
RHZOLD=RHO (I,J) 
RHZ=RHST (I,J) 
ENDIF 
COMPUTE AXIAL VELOCITY, SECOND ESTIMATE 
VZSTT (I, J) = (DW (I, J) * (PIO (I, J-1) -PIO(I, J) ) / (RHZ*PI 
+ * (R(2*I-1)**2-R(2*I-*1)**2))) *VZ(I,J) *RHZOLD/RHZ 
CONTINUE 
CONTINUE 
VELOCITY AT THE BOUNDARY, J-1 
VZSTT(1,1)-DW(1,1)* (-PIO(1, 1) )/ (RHOBI*PI* (R(1) **2-R(3) **2)) 
+ +VZ(1,1) 
VZSTT(M, 1)=DW(M, 1) * (-PIO(M, 1) ) / (RHO (M, 1) *PI*(R(2*M-1)**2)) 
+ +VZ (M, 1) * (RHO(M, 1) /RHST (M, 1)) 
DO 56 I=2,M-1 
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VZ (1,1)=0 
CONTINUE 
RADIAL VELOCITY 
GENERAL FORM OF THE EQUATION 
VRSTT(I,J)- (DS(I,J) * (PIO(I-1,J) -PIO(I, J) *R(I,J) *RHZOLD*AREA) / 
(RHZ*AREA) 
DO 64 I=2,M 
DO 60 J=1,N 
UPWIND TO DETERMINE DENSITY AT VR 
ТЕ (VR(I,J) GD. 0.) THEN 
RHZOLD=RHO (I-1,J) 
RHZ=RHST (I-1, 9) 
ELSE 
RHZOLD=RHO (I, J) 
RHZ=RHST (I,J) 
ENDIF 
COMPUTE RADIAL VELOCITY, SECOND ESTIMATE 
VRSTT (I,J) =(DS(I, J) * (PIO(I-1, J) -PIO(I,J)) /(RHZ*2*PI* 
+ R(2*I-1)*DZ))+VR(I,J) *RHZOLD/RHZ 
CONTINUE 
CONTINUE 
CHECK OVERALL SYSTEM MASS BALANCE 
GI=RHOBI*VZSTT (1, 1) *PI* (R(1) **2-R (3) **2) 
GD=RHST (M, 1) *VZSTT(M, 1) *PI* (R(2*M-1) **2) 
WRITE(3,*) GI,GD 
WRITE(6,*) GI,GD 
RETURN 
END 


SUBROUTINE TO PERFORM THE ENTHALPY BALANCE 


SUBROUTINE ENTHLP (VZSTT, VRSTT, RHO, RHOO, RHST, PST, PIO,H, HO, 
+  UBAR,TSST,TGO, ANN, HST, TGST, AV, QIN, TSSTT) 
COMMON/GEOTRY/R (45) ,EP (20) ,DZ, THKCF (30) , DV (20), DP, DPCF, M,N, DT 
COMMON /BOUCON/PBI, PBO, TBI, RHOBI, HBI 
COMMON / SOLVER/ SOLH (100, 100) 
INTEGER M,N,I,J,K,L,UNITS,NP 
REAL VZSTT(10,11),VRSTT (11,10), RHST (10, 10) , RHO(10, 10) 
REAL PST(10,10),MWH(10,10),MEH(10,10) , MSH(10, 10) , MIH (10, 10) 
REAL AWH(10,10),AEH(10,10),ASH(10,10) , ANH(10, 10) , AOH(10, 10) 
FEET O(10 I0) -PIO(IO,10).H(10,10),HST(10,10), TGST(10, 10) 
REAL CONVH (100) , GAM, HSP, QUAL (10, 10) 
ВЕАТ ВНОО (10, 10) , НО (10,10) 
PEAL ОВАВ (10,10), Т$5Т (10,10), АМ (10), отм (10, 10), ТСО (10, 10) 
REAL ANN(10,10),MIKE(10,10),CONST (10, 10) , TSSTT (10, 10) 
PARAMETER (PI=3.14159, UNITS=1) 
INITIALIZE THE COEFFICIENT MATRIX 
DO 10 I=1,M*N 
DO 5 J=1,N*M 
SOLH (I,J) =0. 
CONTINUE 
CONTINUE 
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CALCULATE ENTHALPY FLUX AND ENTHALPY COEFFICIENTS 
EAST AND WEST FACES 
DO 30 J=2,N 
DO 20 I=1,M 
UPWIND DENSITY 
IF (VZSTT(I,J) .GT. 0.) THEN 
САМ-0.5 
ELSE 
GAM=-0.5 
ENDIF 
WEST FACE MASS FLUX 
MWH (I,J) =((0.5+GAM) *RHST (I, J-1)+(0.5-GAM) *RHST (I,J) ) *PI 
+ * (R(2*I-1) **2-R(2*I+1) **2) *VZSTT (I,J) 
IF (I .EQ. (M-1)) MWH(I,J)=0. 
WEST FACE COEFFICIENT 
AWH (I, J) =MWH (I,J) * (0.5+GAM) 
EAST FACE MASS FLUX 
MEH (I, J-1) =MWH (I,J) 
EAST FACE COEFFICIENT 
AEH (I, J-1) =-MEH(I, J-1) * (0.5-GAM) 
CONTINUE 
CONTINUE 
BOUNDARIES, EAST AND WEST FACES NOT CO-LOCATED 
INLET PORT ASSUMING VELOCITY IS POEITIVE 
MWH (1,1) =RHOBI*VZSTT (1,1) *PI* (R(1)**2-R (3) **2) 
AWH (1,1) =MWH(1, 1) 
OUTLET PORT ASSUMING VELOCITY IS NEGATIVE 
MWH (M, 1) =RHST (M, 1) *VZSTT (M, 1) *PI* (R(2*M-1)**2) 
AWH (M, 1)=0. 
DO 40 I=2,M-1 
MWH (I, 1) =0. 
AWH (I, 1) =0. 
MEH (I,N) =0. 
AEH (I,N) =0. 
CONTINUE 
MEH (1,N) =0. 
AEH (1,N)=0. 
MEH (M,N)=0. 
AEH (M,N)=0. 
NORTH AND SOUTH FACES 
DO 60 J=1,N 
BOUNDARY FACES 
MNH (M, J)=0. 
ANH (M,J)=0. 
MSH (1,J)=0. 
ASH(1,J)=0. 
DO 50 I=2,M 
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UPWINDING DENSITY 
ТЕ (VRSTT(1,J) (СЕ. 0.) THEN 
ALF=0.5 
ELSE 
ALF=-0.5 
ENDIF 
SOUTH FACE MASS FLUX 
MSH (I,J) =((0.5+ALF) *RHST(I-1,J)+(0.5-ALF) *RHST(I,J)) 
+ *VRSTT (I,J) *2*PI*R(2*I-1) *DZ 
SOUTH COEFFICIENT 
ASH(1,J)=MSH(1,J)*(0.5+ALF) 
NORTH FACE MASS FLUX 
MNH (1-1, J) =MSH (I, J) 
NORTH COEFFICIENT 
ANH (1-1, J) =-MNH (1-1, J) * (0.5-ALF) 
CONTINUE 
CONTINUE 
CALCULATE CENTRAL COEFFICIENT AND HEAT INPUT FROM THE SOLID PHASE 
DO 80 J=1,N 
AOH(1,J)= -(AEH(1,J)+AWH(1,J)+ANH(1,J)+ASH(1, 3) +MEH(1,J) 
+ -MWH (1, J) +MNH (1, J) -MSH (1, J)) 
MIKE (1,J)2 (DV(1) *RHST(1, J) /DT-AOH(1, J) ) 
CONST (1,J) 2 (DV(1)/DT) * (PIO(1, J) *RHO(1, J) *HO(1,J)) 
АОН (М, J) = -(АЕН(М, J) +AWH (M, J) +ANH (M, J) +ASH (M, J) +MEH (M, J) 
+ -MWH (M, J) *MNH (M, J) -MSH (M, J} ) 
MIKE (M, J) - (DV (M) *RHST (M, J) /DT-AOH (M, J) ) 
CONST (M, J) = (DV (M) /DT) * (PIO (M, J) *RHO (M, J) *HO (M, J) ) 
DO 70 I-2,M-1 
CPEE-PTCP (PST(I,J), TGO(I,J), UNITS) 


AOH (I, J)= -(AEH(I, J) +AWH (1, J) +ANH (I, J) +ASH (I1, J) +MEH(I, J) 
+ -MWH (I, J) *MNH (I, J) -MSH(I,J)) 

MIKE (I,J)=(DV(I) *RHST(I, J) /DT-AOH (I, J) +UBAR(I, J) *DV(I) 
+ *AV (I) / (CPEE*EP (I) )-((UBAR(I,J) *DV(I) *AV(I) /EP (I) ) **2) 
+ / (ANN (I, J) *CPEE) ) 

CONST (I,J) =UBAR(I, J) * (DV(I) /EP (I) ) *AV(I) * (TSST (I, J) 
+ -TGO (I, J) + (HO(I,J) /CPEE) ) - (C(UBAR(I,J) *DV(I)/EP(I))* 
* AV (I)) **2) * ( (HO(I, J) /CPEE) ) /ANN(I, J) * (DV (I) /DT) 
+ *(PIOGCD,J)TRHO(I,J)*HO(I,J)) 

CONTINUE 
CONTINUE 


CONSTRUCT SOLUTION MATRIX AND CONSTANT VECTOR 
DO 100 I=1,M 
DO 90 J=1,N 
K-2(I-1)*NtJ 
SOLH (K, K) 2-MIKE (I, J) 
ТЕ (J .LI. М) SOLH(K,KTl)*AEH(I,J) 
ТЕ (K+N .LE. M*N) SOLH(K, K+N)=ANK (I, J) 
IE (I .GT. 1) SOLH(K,K-1)-AWH(I,J) 
IF (K-N .GT. 0) SOLH(K,K-N)=ASH(I,J) 
CONVH (K)=-CONST (I,J) 
CONTINUE 
CONTINUE 
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C CALL SOLUTION SUBROUTINES 
C NOTE THAT THE UPDATED ENTHALPIES ARE RETURNED AS A MODIFIED 
C CONVH VECTOR 


NP=100 
CALL LUDCMP (SOLH, (M*N), NP, INDX, D) 
CALL LUBKSB (SOLH, (M*N),NP, INDX, CONVH) 
C CONVERT UPDATED ENTHALPIES TO I,J SUBSCRIPTED FORM 
I-1 
J=1 
DO 101 K=1,M*N 
HST (I, J) =CONVH (K) 
IF ((K-(I-1)*N) .EQ. N) THEN 
J=1 
I-I41 
ELSE 
J=J+1 
ENDIF 
101 CONTINUE 
C UPDATE GAS PHASE TEMPERATURE MATRIX 
DO 103 I=1,M 
DO 102 J=1,N 
TGST(I,J)2PHTEMP (PST(I,J),HST(I, J) , QUAL(I,J) , UNITS) 
102 CONTINUE 
103 CONTINUE 
C UPDATE TSST TO TSSTT AND CALCULATE HEAT TRANSFERRED FROM THE SOLID 
С “TO THE GAS PHASE 
DO 105 I=2,M-1 
DO 104 J=1,N 
TSSTT (I, J) =TSST (I, J) +UBAR (I, J) *DV(I) *AV(I) * (TGST (I,J) 


+ -TGO (I,J) ) / (EP (I) *ANN (I, J) ) 
QIN (I, J) =UBAR (I, J) * (DV (I) /EP (1)) *AV(1I) * (TSSTT (I, J) 
+ =(GST(1,J)) 
104 CONTINUE 
105 CONTINUE 
RETURN 
END 


© 
C THE SECOND MOMENTUM CORRECTOR - THE PRESSURE EQUATION 
с 
SUBROUTINE PRESII (VRST, VZST, VRSTT, VZSTT, P, PST, PO, TGST, RHOO, RHST, 
+  VRSTTT,VZSTTT, PSTT, RHSTT, BVZ, BVR, ARF,AZF,AWR, AER, ASR, ANR, AWZ, 
+ AEZ,ASZ,ANZ, RHO) 
COMMON/GEOTRY/R(45),EP(20),DZ, THKCF (30), DV(20) , DP, DPCF, M, N, DT 
COMMON/BOUCON/PBI,PBO, TBI,RHOBI,HBI 
COMMON/ SOLVER/ SOLC (100, 100) 
INTEGER M,N,I,J,K,UNITS,NP 
REAL VRST(11,10),VRSTT(11,10),V2ST(10,11),VZSTT (10, 11) 
REAL RHOO(10,10),RHST(10,10),P(10,10),PO(10,10), TGST(10, 10) 
REAL HPRI(10,10),HPRIZE(10,10),HPRIZW(10,10),HPRIRN(10,10) 
ВЕАТ НРВТВ$ (10,10), ВНО (10,10) 
КЕАІ, ВУВ (11,10), ВУ (10,11), АВЕ (11,10), АЕ (10,11) 
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REAL DTC(10,10),DTW(10,10),DTE(10,10),DTS(10,10),DTN(10,10) 
REAL VRSTTT (11,10), VZSTTT(10, 11), PIST(10, 10) , RHSTT(10, 10) 
REAL PSTT(10,10),EPS, RHIN, FAEW(10), FAS(10),FAN(10) 
ВЕАГ АУЕГ (10,10), АУЕБИ (10,10), АУЕБЕ (10,10), АУЕ!,$ (10,10) 
REAL AVELN(10,10),CONVC(100),PST(10,10) 
REAL ANR(11,10),ASR(11,10),AER(11,10),AWR(11, 10) 
ВЕДЕ АМА (ТО 11) ASZ(010 11) AB2(10. 11) ЛИТО. ТІ) 
PARAMETER (PI-3.14159, UNITS-1) 
INITIALIZE THE COEFFICIENT MATRIX 
DO 4 I=1,M*N 
DO 2 J=1,M*N 
SOLC (1 J) Ô. 
CONTINUE 
CONTINUE 
GENERAL FORM OF THE EQUATION 
- (DTO-RHST*DV/ (PST*DT))*PIST + DTW*PISTW + DTE*PISTE + DTS*PISTS 
+ DTN*PISTN = ((DV/DT-AO/RHST)**(-1)) * (DEL(H' (VXSTT-VXST) ) *À 
+ DEL (AO ( (RHST-RHO) /RHST) *VXSTT*A) ) -DV/DT* (RHST/PST-RHO/P) 


CALCULATE THE PRESSURE INCREMENT COEFFICIENTS 
EAST AND WEST FACES 
BOUNDARIES 
DO 10 I=2,M-1 
FAEW (I) =PI* (R(2*I-1) **2-R(2*I+1) **2) 
DIW(T,1)0. 
DTW(I,1)=(FAEW(I)**2)*EP(I)/(BVZ(I,1)-(AZF(I,1)/RHST(I,1))) 
CONTINUE 
INLET, ASSUMING THE VELOCITY IS POSITIVE 
FAEW(1)=((R(1)**2-R(3) **2) *PI) 
DTW(1,1)=(FAEW(1) **2) / (BVZ(1,1)-(AZF (1,1) /RHOBI) ) 
OUTLET, ASSUMING THE VELOCITY IS NEGATIVE 
FAEW (M) =((R(2*M-1) **2) *PI) 
DTW (M, 1) x (FAEW (M) **2) / (BVZ (M, 1) - (AZF (M, 1) /RHST (M, 1) )) 
EAST BOUNDARY 
DO 30 1=1,M 
DIE (T,N)=0: 
CO-LOCATED EAST AND WEST FACES 
DO 20 J=2,N 
UPWIND FOR DENSITY 
IF (VZSTT(I,J) .GE. 0.) THEN 
GAM=0 .5 
ELSE 
GAM=-0.5 
ENDIF 
WEST AND EAST PRESSURE INCREMENT COEFFICIENT 
RHIN= (0.5+GAM) *RHST (I,J-1)+(0.5-GAM) *RHST (I, J) 
DIW (I, J) = (FAEW (I) **2) *EP (I) / (BVZ (I, J) - (AZF (I, J) /RHIN) ) 
DTE (I, J-1) =DTW(I,J) 
CONTINUE 
CONTINUE 
NORTH AND SOUTH FACE COEFFICIENTS 
BOUNDARY CONDITIONS 
CALCULATE THE FACE AREAS 
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DO 50 I=2,M 
FAS(I)=(2*PI*R(2*I-1)*DZ) 
FAN (I-1) =FAS (I) 


CONTINUE 
FAN (M) =0. 
DO 70 J=1,N 
DTN (M, J) =0. 
DTS(1,J)=0. 


CO-LOCATED NORTH AND SOUTH FACE COEFFICIENTS 
DO 60 I=2,M 
UPWIND DENSITY 
IF (VRSTT(I,J) .GE. 0.) THEN 
АТЕ=О 5 
ELSE 
ALF=-0.5 
ENDIF 
EPS=MIN (EP (I),EP(I-1)) 
NORTH AND SOUTH COEFFICIENTS 
RHIN=(0.5+ALF) *RHST (I-1,J)+(0.5-ALF) *RHST (I,J) 
DTS (I,J) =(FAS (I) **2) *EPS/ (BVR(I,J)- (ARF (I, J) /RHIN) ) 
DTN (I-1,J)=DTS (I, J) 
CONTINUE 
CONTINUE 
CENTRAL POINT COEFFICIENT 
DO 90 I=1,M 
DO 80 J=1,N 
DTC (I1, J) =DTE (I, J) +DIW (I, J) +DTS (I, J) +DTN (I,J) 
+ + ( (RHST (I, J) *DV(I))/ (PST(I,J) *DT)) 
CONSTRUCT THE COEFFICIENT MATRIX 
K= (1-1) *N+J 
SOLC(K,K)s-DTC(I,J) 
IF (J .LT. N) SOLC(K, K+1)=DTE(I,J) 
IF ((K+N) .LE. (M*N)) SOLC(K, K+N)=DTN (I,J) 
IF (J .GT. 1) SOLC(K, K-1) =DTW(I,J) 
IF ((K-N) .GT. 0) SOLC(K, K-N)=DTS (I,J) 
CONTINUE 
CONTINUE 
THE CONSTANT TERMS 
DO 110 I=2,M-1 
DO 100 J=2,N-1 
UPWIND FOR DENSITY 
IF (VZSTT(1,J+1) .GE. 0.) THEN 
RHE=RHST (I,J) 
RHEE- (RHST (I,J) -RHO(I, J) )/RHO(I, J) 
ELSE 
RHE=RHST (I, J+1) 
RHEE= (RHST (I, J+1) -RHO(I, J+1) ) /ВНО(Т, 9+1) 
ENDIF 
ТЕ (VRSTT(I+1,J) .GE. 0.) THEN 
RHN=RHST (I,J) 
RHNN= (RHST (I,J) -RHO(I,J)) /RHO(I,J) 
ELSE 
RHN=RHST (I+1,J) 
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RHNN= (RHST (I+1, J) -RHO(I+1, J) ) /RHO(I+1, J) 
ENDIF 
Е CLACULATION OF THE CONSTANTS ON CO-LOCATED FACES 
HPRIZE(I,J)=(AWZ (I, J+1) * (VZSTT(I, J) -VZST(I, J) ) 
+AEZ (1,3+1)* (VZSTT(1,J+2)-VZST(1,J+2)) 
+ANZ (I, J+1) * (VZSTT (I4+1, J+1) -VZST (I4+1, J+1) ) 
+ASZ (I, J+1)* (VZSTT(I-1,J+1)-VZST(I-1,J+1))) 
*FABW(I) / (BVZ (I, J+1) -(AZF (I, J+1) /RHE) ) 
HPRIZW(I,J+1) =HPRIZE (I,J) 
HPRIRN(I,J)=(AWR(I+1, J) * (VRSTT (I+1, J-1) -VRST(I+1, J-1) ) 
+AER (14+1,J3)* (VRSTT(1+1,J+1)-VRST(1+1,J+1)) 
+ASR (I+1, J3) * (VRSTT (I, J3) -VRST (I, 3)) 
+ANR (I+1,J) * (VRSTT (I+2, J) -VRST (I+2,J))) 
*FAN (I) / (BVR(I+1,J) - (ARF (I+1,J)/RHN) ) 
HPRIRS(I+1, J) =HPRIRN(I, J) 
AVELE (I, J)=AZF (I, J+1) *RHEE*VZSTT (I, J+1) 
+ *FAEW (I) / (BVZ(I, J+1) -(AZF (I, J+1) /RHE) ) 
AVELW (I, J+1)=AVELE (I, J) 
AVELN (I, J) =ARF (I+1, J) *RHNN* VRSTT (I+1, J) 


+ + + + 


+ + + + 


+ «ҒАМ (І) / (BVR(I+1, J) - (ARF (I+1,J) /RHN) ) 
AVELS (I+1,J)=AVELN (I, J) 
100 CONTINUE 
110 CONTINUE 


C CALCULATION OF THE CONSTANTS ALONG THE BOUNDARIES 
C THE INLET PORT 
TE (М26ТТ(1,2) СЕ 0.) THEN 
RHE=RHST (1,1) 
RHEE- (RHST (1, 1) -RHO (1, 1) )/RHO(1, 1) 
ELSE 
RHE=RHST (1, 2) 
RHEE- (RHST (1,2) -RHO(1,2) )/ RHO(1, 2) 
ENDIF 
HPRIZE (1,1)=(AWZ (1,2) * (VZSTT(1,1)-VZST(1,1)) 
+AEZ (1,2) *(VZSTT(1,3) -VZST(1, 3)) 
+ANZ (1,2) * (VZSTT(2,2)-VZST(2,2))) 


+ *FAEW(1)/(BVZ(1,2)- (AZF (1,2) / RHE)) 
HPRIZW(1,1)-2(AEZ(1,1) *(VZSTT (1,2) -VZST(1,2))) 
+ *FAEW (1) / (BVZ (1,1) - (AZF (1, 1) /RHOBI)) 


IF (VRSTT (2,1) .GE. 0.) THEN 
RHN=RHST (1,1) 
RHNN- (RHST (1, 1) -RHO(1,1)) /RHO(1, 1) 
ELSE 
RHN-RHST (2,1) 
RHNN- (RHST (2, 1) -RHO (2,1) ) /RHO (2, 1) 


ENDIF 

HPRIRN(1,1)-(AER(2,1) * (VRSTT (2,2) -VRST(2,2)) 
+ +ANR (2, 1) * (VRSTT (3,1) -VRST(3,1))) 
*FAN (1) / (BVR(2,1)-(ARF(2,1)/RHN)) 


HPRIRS (1,1)=0. 
AVELE (1,1)=AZF (1, 2) *RHEE*VZSTT (1, 2) 
+ *FAEW (1) / (BVZ (1, 2) = (AZF (1,2) /RHE) ) 
AVELW (1,1) =0. 
AVELN (1,1) =ARF (2, 1) *RHNN*VRSTT (2, 1) 
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*FAN (1) / (BVR(2, 1) - (ARF (2, 1) / RHN) ) 
AVELS (1,1)=0. 


С SHE OUTLET PORT 


+ 


+ 


ТЕ OROUASTT(M 2)..GE. 0.) THEN 
ВНЕ=ВН$Т (М, 1) 
RHEE- (RHST (M, 1) -RHO (M, 1)) /RHO (M, 1) 
ELSE 
RHE=RHST (M, 2) 
RHEE= (RHST (M, 2) -RHO (M, 2) ) /RHO (M, 2) 
ENDIF 
IF (VRSTT(M,1) .GE. 0.) THEN 
RHS=RHST (M-1,1) 
RHSS- (RHST (M-1, 1) -RHO (M-1,1)) /RHO (M-1,1) 
ELSE 
RHS=RHST (M, 1) 
RHSS= (RHST (M, 1) -RHO (M, 1) ) /RHO (M, 1) 
ENDIF 
HPRIZE (M, 1) =(AWZ (M, 2) * (VZSTT (M, 1) -VZST (M, 1)) 
*AEZ (M, 2) * (VZSTT (M, 3) -VZST (M, 3)) 
+ASZ(M, 2) * (VZSTT (M-1, 2) -VZST(M-1, 2) )) 
*FAEW (M) / (BVZ (M, 2) - (AZF (M, 2) / RHE)) 
HPRIZW(M, 1)- (AEZ (M, 1) * (VZSTT (M, 2) -VZST (M, 2) )) 
*FAEW (M) / (BVZ (M, 1) - (AZF (M, 1) /RHST (M, 1))) 
HPRIRN (M, 1)=0. 
HPRIRS (M, 1) - (AER (M, 1) * (VRSTT (M, 2) -VRST (M, 2)) 
+ASR (M, 1) * (VRSTT (M-1, 1) -VRST (M-1, 1))) 
*FAS (M) / (BVR (M, 1) - (ARF (M, 1) /RHS)) 
AVELE (M, 1) -AZF (M, 2) *RHEE*VZSTT (M, 2) 
*FAEW (M) / (BVZ (M, 2) - (AZF (M, 2) / RHE)) 
AVELW (M, 1) -AZF (M, 1) * ( (RHST (M, 1) -RHO (M, 1) ) /RHO (M, 1) ) *VZSTT (M, 1) 
*FAEW (M) / (BVZ (M, 1) - (AZF (M, 1) /RHST (M, 1))) 
AVELN (M, 1) =0. 
AVELS (M, 1) -ARF (M, 1) *RHSS*VRSTT (M, 1) 
*FAS (M) / (BVR (M, 1) - (ARF (M, 1) /RHS)) 


E THE CORNER CONTROL VOLUMES 


IF (VZSTT(1,N) .GE. 0.) THEN 
RHW=RHST (1,N-1) 
RHWW= (RHST (1,N-1)-RHO(1,N-1))/RHO(1,N-1) 
ELSE 
RHW=RHST (1,N) 
RHWWz (RHST (1,N) -RHO(1,N) )/RHO(1,N) 
ENDIF 
IF (VRSTT(2,N) .GE. 0.) THEN 
RHN=RHST (1,N) 
RHNN- (RHST (1,N) -RHO(1,N) ) /RHO(1, N) 
ELSE 
RHN-RHST (2, N) 
RHNNs (RHST (2, N) -RHO(2,N) ) /RHO (2, N) 
ENDIF 
HPRIZE(1,N)=0. 
HPRIZW(1,N) =(AWZ(1,N) * (VZSTT(1,N-1)-VZST(1,N-1)) 
+ANZ(1,N) *(VZSTT(2,N) -VZST(2,N))) 
*FAEW (1)/ (BVZ (1, N) - (AZF (1, N) / RHW) ) 
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HPRIRN(1,N)-(AWR(2,N)* (VRSTT (2, N-1) -VRST (2, N-1)) 
+ +ANR(2,N) * (VRSTT(3,N) -VRST(3,N))) 
+ *FAN(1)/(BVR(2,N)-(ARF(2,N)/RHN)) 
HPRIRS(1,N)=0. 
AVELE (1,N)=0. 
AVELW (1,N)=AZF (1,N) *RHWW*VZSTT(1,N) 


+ *FAEW(1)/(BVZ(1,N)- (AZF (1, N) / RHW) ) 
AVELN (1,N)-ARF (2,N) *RHNN*VRSTT (2, N) 
+ *FAN(1)/(BVR(2,N)-(ARF(2,N)/RHN)) 


AVELS (1,N)=0. 
IF (VZSTT(M,N) .GE. 0.) THEN 
RHW=RHST (M, N-1) 
RHWW- (RHST (M, N-1) - RHO (M, N21) ) /RHO (M, N-1) 
ELSE 
RHW=RHST (M, N) 
RHWW- (RHST (M, N) -RHO (M, N) ) /RHO (M, N) 
ENDIF 
IF (VRSTT(M,N) .GE. 0.) THEN 
RHS=RHST (M-1,N) 
RHSS- (RHST (M-1,N) -RHO (M-1,N)) /RHO (M-1,N) 
ELSE 
ВН$=ВН$Т (М, М) 
RHSS- (RHST (M, N) -RHO (M, N) )/RHO (M, N) 
ENDIF 
HPRIZE (M,N) =0. 
HPRIZW(M, N) = (AWZ (M,N) * (VZSTT (M, N-1) -VZST (M, N-1) ) 
+ +ASZ (M,N) * (VZSTT (M-1,N) -VZST (M-1,N))) 
+ *FAEW (M) / (BVZ (M, N) - (AZF (M, N) / RHW) ) 
HPRIRN (M,N) =0. 
HPRIRS (M, N) = (AWR (M, N) * (VRSTT (M, N-1) -VRST(M,N-1)) 
Б +ASR(M, N) * (VRSTT (M-1,N)-VRST(M-1,N))) 
+ *FAS (M) / (BVR (M, N) - (ARF (M, N) /RHS)) 
AVELE (M,N)=0. 
AVELW (M, N) =AZF (M, N) *RHWW*VZSTT (M, N) 
+ *FAEW (M) / (BVZ (M, N) - (AZF (M, N) / RHW) ) 
AVELN (M,N)=0. 
AVELS (M, N) =ARF (M, N) *RHSS*VRSTT (M,N) 
+ *FAS (M) / (BVR(M,N) - (ARF (M, N) /RHS) ) 
E THE EAST AND WEST BOUNDARY CONTROL VOLUMES 
DO 120 I=2,M-1 
ТЕ (VZSTT(1,2) .GE. 0.) THEN 
RHE=RHST (1,1) 
RHEE- (RHST(I,1) -RHO(I,1)) /RHO(I,1) 
ELSE 
RHE-RHST(I, 2) 
RHEE- (RHST(I,2) -RHO(I,2)) /RHO (I,2) 
ENDIF 
IF (VRSTT(I+1,1) .GE. 0.) THEN 
RHN=RHST (I, 1) 
RHNN- (RHST(I, 1) -RHO(I,1))/RHO(I,1) 
ELSE 
RHN=RHST (I+1,1) 
RHNN=(RHST(I+1,1)-RHO(I+1,1))/RHO(I+1,1) 
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+ 


ENDIF 

IF (VRSTT(I,1) .GE. 0.) THEN 
RHS=RHST(I-1, 1) 
RHSS- (RHST(I-1,1) -RHO(I-1,1)) /RHO(I-1,1) 
ELSE 
RHS-RHST(I, 1) 
RHSS- (RHST (I, 1) -RHO(I,1)) /RHO(I, 1) 

ENDIF 

HPRIZE(I,1)-(AEZ(I,2)*(VZSTT(I,3) -VZST(I,3)) 
-ANZ(I,2)*(VZSTT(I-*1,2)-VZST(I-*1,2)) 
+ASZ (1,2) * (VZSTT(I-1,2)-VZST(I-1,2))) 
*FAEW(I)/(BVZ(I,2) -(AZF(I,2)/RHE)) 

HPRIZW(I,1)=0. 

HPRIRN (I, 1)=(AER(I+1,1) * (VRSTT (I+1,2) -VRST (I+1, 2) ) 
+ANR (I+1,1) * (VRSTT (I+2,1)-VRST (I+2,1)) 
+ASR(I+1,1)* (VRSTT(I,1)-VRST(I,1))) 
*FAN (I) / (BVR(I+1,1)- (ARF (I+1,1)/RHN) ) 

HPRIRS (I,1)=(AER(I,1) * (VRSTT(I,2) -VRST(I,2) ) 
+АМВ (Т,1) * (VRSTT (1+1,1)-VRST(1+1,1)) 
+ASR (1,1) * (VRSTT (1-1,1)-VRST(I-1,1))) 
*FAS (1) / (BVR(1,1)-(ARF (1,1)/RHS)) 

AVELE (1,1)=AZF (1,2) *RHEE*VZSTT (1,2) 
*FAEW(I)/(BVZ(I,2)- (AZF(I,2) /RHE)) 

AVELW(I,1)-0. 

AVELN(I,1)-ARF(I-1,1)*RHNN*VRSTT (I-*1,1) 

*FAN (I)/ (BVR(I*1,1)- (ARF (I*1,1)/RHN)) 

AVELS(I,1)-ARF(I,1)*RHSS*VRSTT (I, 1) 

*FAS (I) / (BVR(I,1)- (ARF (I, 1) /RHS) ) 

IF (VZSTT(I,N) .СЕ. 0.) THEN 

RHW=RHST (I, N-1) 

RHWW- (RHST (I,N-1) -RHO (I,N-1)) /RHO(I,N-1) 
ELSE 

RHW-RHST (I,N) 

RHWW- (RHST (I,N) -RHO(I,N))/RHO(I,N) 

ENDIF 

IF (VRSTT(I,N) .GE. 0.) THEN 

RHS=RHST (I-1,N) 

RHSS- (RHST(I-1,N) -RHO(I-1,N)) /RHO (I-1,N) 
ELSE 

RHS-RHST(I,N) 

RHSS- (RHST (I, N) -RHO (I,N)) /RHO (I,N) 

ENDIF 

IF (VRSTT(I+1,N) .GE. 0.) THEN 

RHN=RHST (1,N) 

RHNN- (RHST (I, N) -RHO (I,N)) 
ELSE 

RHN=RHST (I+1,N) 

RHNN= (RHST(I+1,N) ~RHO(I+1,N) ) 

ENDIF 
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+ 


+ 


120 


HPRIZE(1,N)=0. 

HPRIZW(I,N)=(AWZ (I,N) * (VZSTT(I,N-1)-VZST(I,N-1)) 
+ANZ (I,N) * (VZSTT (I+1,N) -VZST(I+1,N) ) 
+ASZ (I,N) * (VZSTT (I-1,N) -VZST(I-1,N))) 
*FAEW(I)/(BVZ(I,N)-(AZF(I,N)/RHW)) 

HPRIRN (I,N)=(AWR(I+1,N) * (VRSTT (I+1,N-1) -VRST(I+1,N-1)) 
+ANR (I+1,N) * (VRSTT(I+2,N)-VRST (I+2,N)) 
+ASR(I+1,N) * (VRSTT(I,N)-VRST(I,N))) 

*FAN (1) / (BVR(I+1,N)- (ARF (I+1,N)/RHN) ) 

HPRIRS (I,N) =(AWR(I,N) * (VRSTT(I,N-1)-VRST(I,N-1)) 
+ANR (I,N) * (VRSTT (I+1,N) -VRST(I+1,N) ) 
+ASR(I,N)*(VRSTT(I-1,N)-VRST(I-1,N))) 
*FAS(I)/(BVR(I,N)-(ARF(I,N)/RHS)) 

АУЕҺЕ(І,М)-0. 

AVELW(I,N)-AZF(I,N)*RHWW*VZSTT (I,N) 
*FAEW(I)/(BVZ(I,N)-(AZF(I,N) /RHW)) 

AVELN(I,N)-ARF(I-1,N)*RHNN*VRSTT (I41,N) 

*FAN (I)/ (BVR(I+1,N)- (ARF (I+1,1)/RHN) ) 

AVELS(I,N)=ARF (I,N) *RHSS*VRSTT (I,N) 
*FAS(I)/(BVR(I,N)-(ARF(I,N)/RHS)) 


CONTINUE 


С NORTH AND SOUTH BOUNDARY CONTROL VOLUMES 


- 


+ 


DO 130 J=2,N-1 


IF (VZSTT(M,J) .GE. 0.) THEN 
RHW-RHST (M, J-1) 
RHWW» (RHST (M, J-1) -RHO (M, J-1) ) /RHO (M, J-1) 
ELSE 
RHW=RHST (M, J) 
RHWW= (RHST (M, J) -RHO (M, J) ) /RHO (M, J) 
ENDIF 
IF (VZSTT(M,J+1) .GE. 0.) THEN 
RHE=RHST (M, J) 
RHEE- (RHST (M, J) -RHO (M, J) ) /RHO (M, J) 
ELSE 
RHE-RHST (M, J*1) 
RHEE-s (RHST (M, J* 1) -RHO (M, J*1) )/RHO (M, J*1) 
ENDIF 
IF (VRSTT(M,J) .GE. O.) THEN 
RHS=RHST (M-1, J) 
RHSS- (RHST (M-1, J) - RHO (M-1,J)) /RHO (M-1,J) 
ELSE 
RHS-RHST (M, J) 
RHSS- (RHST (M, J) - RHO (M, J) ) /RHO (M, J) 
ENDIF 
HPRIZE(M, J) = (AWZ (M, J+1) * (VZSTT (M, J) -VZST (M, J) ) 
+AEZ (M, J+1) * (VZSTT (M, J+2) -VZST (M, J+2) ) 
+ASZ (M, J+1) * (VZSTT (M-1,J+1) -VZST(M-1,J+1))) 
*FAEW(M) / (BVZ (M, J+1) - (AZF (M, J+1) /RHE) ) 
HPRIZW(M,J) = (AWZ (M, J) * (VZSTT (M, J-1) -VZST (M,J-1)) 
+AEZ (M, J) * (VZSTT (M, J+1) -VZST (M, J+1) ) 
+ASZ (M, J) * (VZSTT (M-1,J) -VZST (M-1,J))) 
*FAEW (M) / (BVZ (M, J) - (AZF (M, J) /RHW)) 
HPRIRN(M,J)=0. 
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AVELE (M, J) =AZF (M, J+1) *RHEE*VZSTT (M, J+1) 


+ *FAEW (M) / (BVZ (M, J* 1) - (AZF (M, J*1) /RHE) ) 
AVELW (M, J) =AZF (M, J) *RHWW*VZSTT (M, J) 
+ *FAEW (M) / (BVZ (M, J) - (AZF (M, J) /RHW)) 


АУЕГМ (М, 9) =0. 
ТЕ (VZSTT(1 J) GE. 0.) “THEN 
RHW=RHST(1,J-1) 
RHWW- (RHST (1, J-1) -RHO(1,J-1)) /RHO(1,J-1) 
ELSE 
RHW=RHST (1, J) 
RHWW= (RHST (1, J) -RHO(1,J)) /RHO(1,J) 
ENDIF 
IF (VZSTT(1,J+1) .GE. 0.) THEN 
RHE=RHST (1, J) 
RHEE- (RHST (1, J) -RHO(1,J)) /RHO(1,J) 
ELSE 
RHE=RHST (1, J+1) 
RHEE= (RHST(1,J+1) -RHO(1,J*1)) /RHO(1,J41) 
ENDIF 
IF (VRSTT(2,J) .GE. 0.) THEN 
RHN-RHST(1,J) 
RHNNz (RHST (1, J) -RHO (1,J)) /RHO(1,J) 
ELSE 
RHN-RHST (2, J) 
RHNN- (RHST (2, J) -RHO (2, J) ) /RHO(2,J) 


ENDIF 
HPRIZE(1,J)=(AWZ (1, J+1) * (VZSTT (1,3) -VZST (1,J)) 

T +AEZ (1, J0+1) * (VZSTT (1, J+2) -VZST (1, J+2) ) 

Е +ANZ (1,J+1) * (VZSTT (2, J+1) <VZST (2, J+1))) 

+ *FAEW(1)/(BVZ(1,J*1) - (AZF (1, J*1) /RHE)) 
HPRIZW(1,J)=(AWZ (1,3) * (VZSTT (1,J-1)-VZST(1,J-1)) 

$ +AEZ (1,J) * (VZSTT (1, J+1) -VZST (1, J+1)) 

+ +ANZ (1,3) * (VZSTT (2,7) -VZST(2,J))) 

+ *FAEW(1)/ (BVZ (1, J) -(AZF (1, J) /RHW) ) 
HPRIRN(1,J)- (AWR (2, J) * (VRSTT (2,J-1) -VRST (2,J-1)) 

$ +AER (2, J) * (VRSTT (2, J+1) -VRST (2, J+1) ) 

} +ANR (2, J) * (VRSTT (3, J) -VRST(3,J))) 

+ *FAN (1) / (BVR (2, J) - (ARF (2, J) /RHN)) 
HPRIRS(1,J)-20. 
AVELE (1, J) =AZF (1,J+1) *RHEE*VZSTT (1, J4*1) 

+ «РАЕН (1) / (BVZ (1, J+1) -(AZF (1, J+1) /RHE) ) 
AVELW(1,J) =AZF (1,7) *RHWW*VZSTT (1, J) 

+ *PAEW (1) / (BVZ (1, J) - (AZF (1, J) /RHW) ) 
AVELN (1, J) =ARF (2, J) *RHNN*VRSTT (2, J) 

+ *FAN (1) / (BVR (2, J) - (ARF (2, J) /RHN) ) 
AVELS(1,J)-0. 

130 CONTINUE 


С CALCULATE THE GRADIENT OF THE CONSTANT TERMS ACROSS THE CONTROL 
С VOLUME 
DO 150 I=1,M 
DO 140 J=1,N 
HPRI (I,J) =HPRIZE(I,J) -HPRIZW(I, J) +HPRIRN (I,J) -HPRIRS (I,J) 
AVEL (I,J) =AVELE (I,J) -AVELW(I, J) +AVELN (I, J) -AVELS (I, J) 
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140 CONTINUE 
150 CONTINUE 
C CONSTRUCT THE CONSTANT VECTOR 
DO 170 I=1,M 
DO 160 J=1,N 
K-(I-1)*N4J 
CONVC (K) EHPRI (I, J) -AVEL (I,J) * (DV(I) /DT) * 


+ ((RHST(I,J)/PST(I,J))- (RHOO(I, J) /PO(I,J))) 
160 CONTINUE 
1779 CONTINUE 
© CALL SOLUTION ROUTINES 
NP=100 
NR=M*N 


CALL LUDCMP (SOLC,NR,NP, INDX,D) 
CALL LUBKSB(SOLC,NR,NP, INDX, CONVC) 
C | SOLUTION RETURNED AS A MODIFIED CONVC VECTOR 
I-1 
J=1 
DO 180 K=1,NR 
PSTT (I,J) =PST(I, J) +CONVC (K) 
PIST (I,J) =CONVC (K) 
IF ((K-(I-1)*N) .EQ. N) THEN 
J=1 
I=I+1 
ELSE 
J=J+1 
ENDIF 
180 CONTINUE 
C UPDATE DENSITIES 
DO 200 I=1,M 
DO 190 J=1,N 
RHSTT(I,J)=PTDENS (PSTT(I,J),TGST(I,J),UNITS) 
190 CONTINUE 
200 CONTINUE 
UPDATE VELOCITIES 
AXIAL VELOCITY FIELD 
GENERAL FORM OF THE EQUATION 
VZSTTT = (HPRIZW-DTW* (PIST(I,J)-PIST(I,J-1))-AVELW*VZSTT) 
* (1/ (RESTT*FAEW) ) + (RHOLD*VZSTT) /RHZ 


Q QQ Q O Q OQ O 


UPWIND FOR DENSITY 
DO 220 I=1,M 
C IF (I .EQ. M-1) GOTO 220 
DO 210 J=2,N 
IF (VZSTT(I,J) .GE. 0.) THEN 
RHZ=RHSTT (I, J-1) 
RHOLD=RHST (I, J-1) 
ELSE 
RHZ=RHSTT (I,J) 
RHOLD=RHST (I,J) 
ENDIF 
C COMPUTE AXIAL VELOCITY FIELD 
VZSTTT(I,J)=(HPRIZW(I,J)-(DTW(I,J)*(PIST(I,J)-PIST(I,J-1) 
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+ ))= (AVELW(1,J)*V2STT(I,J)))*(1/ (RH2*FAEW(I)))+(RHOLD* 
T VZSTT(I,J))/RHZ 
210 CONTINUE 
220 CONTINUE 
e COMPUTE INLET VELOCITY ASSUMING IT IS POSITIVE 
VZSTTT(1,1)-(HPRIZW(1,1)-(DTW(1,1)*PIST(1,1))-(AVELW(1,1) 
+ *VZSTT(1,1)))* (1/(RHOBI*FAEW (1)))+VZSTT(1,1) 
Ç COMPUTE OUTLET VELOCITY ASSUMING IT IS NEGATIVE 
VZSTTT (M, 1) =(HPRIZW(M,1)-(DTW(M, 1) *PIST (M, 1) ) - (AVELW(M, 1) 
+ *VZSTT (M, 1)))* (1/ (RHSTT (M, 1) *FAEW(M) ))+(RHST(M, 1) *VZSTT(M, 1) 
+ /RHSTT (M, 1) ) 
COMPUTE RADIAL VELOCITY FIELD 
GENERAL FORM OF THE EQUATION 
VRSTTT - (HPRIRS-DTS*(PIST(I,J)-PIST(I-1,J)) -AVELS*VRSTT) 
* (1/ (RHSTT*FAS) ) +RHOLD*VRSTT/RHZ 


O O O OO о 


UPWIND FOR DENSITY 
DO 240 I=2,M 
DO 230 J=1,N 
IF (VRSTT(I,J) .GT. 0.) THEN 
RHOLD=RHST (I-1,J) 
RHZ=RHSTT (I-1, J) 
ELSE 
RHOLD=RHST (I,J) 
RHZ=RHSTT (I,J) 
ENDIF 
C COMPUTE THE RADIAL VELOCITY FIELD 
VRSTTT (I, J) =(HPRIRS (I,J) - (DTS (1, J) * (PIST(I, J) -PIST(I=-1, J) 


+ ) ) ~(AVELS (I, J) *VRSTT (I,J))) * (1/ (RHZ*FAS (I) )) 
+ + (RHOLD* VRSTT (I, J) /RHZ) 

230 CONTINUE 

240 CONTINUE 


С CHECK OVERALL SYSTEM MASS BALANCE 
GI=RHOBI*VZSTTT (1,1) *FAEW (1) 
GD=RHSTT (M, 1) *VZSTTT (M, 1) *FAEW (M) 


с WRITE (3,*) GI,GD 
WRITE (6,*) GI,GD 
RETURN 
END 
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Cm» Cp 


EM 


f(r,t) 


Go 


h 


NOMENCLATURE 


Area of control volume face 

Finite difference coefficients for momentum balance 
Area variable for surface integrals 

Particle surface area per unit volume of the bed 
Specific heat of the gas at bulk temperature 
Specific heat of the gas at the film temperature 
Average fuel particle specific heat 

Finite difference coefficient for pressure equation 
Equivalent pipe diameter 

Fuel particle diameter 

Emissivity 

General transport function 

Friction factor 

Force 

Superficial mass velocity 

Specific enthalpy 

Heat transfer coefficient 

Total enthalpy of a control volume 


Finite difference operator, The sum of the E, N, S, W control 
volume variables multiplied by their respective coefficients 


Colburn j factor 
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ЕШ! n Ш 


) 


Ө 


S 


State Variable 


~ 


530807 53 


Si 


Vo 


М 


Conductivity 
Radiation conductivity across a void space to an adjacent particle 


Effective bed thermal conductivity 

Length of control volume 

Mass 

Mass flow rate 

Fuel particle mass per unit particle surface area 
Mass flow through a control volume face 

Unit normal vector 

Pressure 

Heat generation per unit particle outside surface area 
Heat input per unit time 

Reynolds number 

Surface area 

Temperature, Pressure, Density or Enthalpy 
Time 

‘Temperature 

Effective fuel particle temperature 


Effective overall heat transfer coefficient for solid to gas heat 
transfer 


Superficial velocity component in the direction of the momentum 
balance 


Component of the superficial velocity perpendicular to the face 
of the control volume 


Superficial velocity magnitude 
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Va vz 
и j 
En 


Greek Letters 


u 


Subscripts 
a 


су 


о 


Volume 
Axial velocity 
Radial velocity 


Length variable 


Upwinding factor 

Axial momentum correction factor 

Axial momentum correction factor 

Void fraction 

Overall momentum correction factor 

Ratio of solid phase conductivity to gas phase conductivity 
Viscosity 

Viscosity of the coolant at film temperature 

Density 

Heat transfer time constant 


Dummy variable used when equation applies to more than one 
variable; in heat effective thickness of gas film between particles 


Particle shape factor 


Axial 
Control volume 


Identifiers for the faces of a control volume, align in accordance 
with the compass 


Identifiers for the balance variable in the adjacent control vol- 
umes, align in accordance with the compass 
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E 





С 

Р 

$ 

x 
Superscripts 
n ntl 

x 


Identifier for balance variable in the central control volume 
Gas Phase 

Radial. In heat transfer, radiation 

Solid Phase 


Dummy variable for operations that apply to several different 
subscripts such as the directional face identifiers (e, w, n, s) 


Axial direction 


Time step identifiers 


Estimate of time n+1 value of a variable, number of * defines 
which estimate 
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А-2 


В-1 


В-2 


В-3 


В-4 


B-5 
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